LakeCC: a tool for efficiently identifying lake basins with application to palaeogeographic reconstructions of North America

Along the margins of continental ice sheets, lakes formed in isostatically depressed basins during glacial retreat. Their shorelines and extent are sensitive to the ice margin and the glacial history of the region. Proglacial lakes, in turn, also impact the glacial isostatic adjustment due to loading, and ice dynamics by posing a marine‐like boundary condition at the ice margin. In this study we present a tool that efficiently identifies lake basins and the corresponding maximum water level for a given ice sheet and topography reconstruction. This algorithm, called the LakeCC model, iteratively checks the whole map for a set of increasing water levels and fills isolated basins until they overflow into the ocean. We apply it to the present‐day Great Lakes and the results show good agreement (∼1−4%) with measured lake volume and depth. We then apply it to two topography reconstructions of North America between the Last Glacial Maximum and the present. The model successfully reconstructs glacial lakes such as Lake Agassiz, Lake McConnell and the predecessors of the Great Lakes. LakeCC can be used to judge the quality of ice sheet reconstructions. © 2019 The Authors Journal of Quaternary Science Published by John Wiley & Sons Ltd.


Introduction
During the retreat of ice sheets, vast amounts of freshwater are released. The retreating ice also leaves a deeply depressed topography that only slowly recovers due to delayed glacial isostatic adjustment (GIA). In this depressed terrain along the ice margin, meltwater accumulates and forms lakes. These can be the size of an inland sea and impose marine-like boundary conditions on the ice sheet (Tweed and Carrivick, 2015), which may impact the ice dynamics. The largest proglacial lake in North America after the Last Glacial Maximum (LGM) was Lake Agassiz, which formed along the southern margin of the Laurentide Ice Sheet (LIS) (Teller and Leverington, 2004). In Europe, when the Fennoscandian Ice Sheet retreated, the Baltic Ice Lake formed in the basin of the Baltic Sea (Björck, 1995). Today, strandlines and lake sediments provide evidence of the extent of these ancient lakes.
Due to retreating ice and ongoing GIA, the basins evolved over time. The opening of a new, lower outlet could lead to the sudden drainage of a proglacial lake. These events can be observed today with lakes next to alpine glaciers (Walder et al., 2006), but it also happened with ancient lakes. One example is the catastrophic drainage of Lake Missoula, that, due to repeated advances and retreats of the ice sheet, filled and drained more than 40 times (Waitt, 1985). These outburst floods eroded the Channeled Scablands in Washington (Bretz, 1923;Pardee, 1942;Waitt, 1985;Hanson et al., 2012).
Many studies have shown the impact of freshwater hosing on ocean circulation and global climate (e.g. Manabe and Stouffer, 1997;Lohmann and Schulz, 2000;Kageyama et al., 2009). In this regard the evolution of proglacial lakes is of particular relevance. Many hypotheses about the causes of sudden cooling events invoke the (partial) drainage of glacial lakes into the ocean (e.g. Broecker et al., 1989;Barber et al., 1999;Broecker, 2006;Murton et al., 2010;Carlson and Clark, 2012;Wickert, 2016).
The load of the lakes' mass also impacts GIA and should therefore ideally also be considered when reconstructing palaeo ice sheets, as done by Tarasov and Peltier (2006) and Lambeck et al. (2010). Although the impact on GIA is relatively small compared to the ice load, this secondary effect will improve the results of palaeogeography reconstructions. By comparing reconstructed lakes with geological data from the field, these can be fine-tuned.
In this study we present the LakeCC model (CC stands for 'connected components') that locates closed basins that would fill with water for given ice and terrain maps. It is applied to reconstructions of the North American ice-complex after the LGM to the present day (PD) from the NAICE (Gowan et al., 2016b) and ICE-6G (Peltier et al., 2015) models. The results of the latter, however, are only briefly discussed here. After a description of the LakeCC model and its application, the results are compared with geological data. Finally, we discuss the results with regard to possible improvements of the ice margins and different drainage scenarios of Lake Agassiz.

Bed topography
As with most hydrological applications where water is routed through a digital elevation model (DEM), spatial resolution is a critical factor. Narrow spillways can have a significant impact on the results. Increasing resolution might add accuracy but also brings big challenges with regard to computational resources and suitable numerical implementation that is capable of efficiently processing big data (Barnes, 2016).
Although high resolution might be a critical issue, it might not always be possible to obtain a data set of sufficient quality, especially when using reconstructed palaeo data sets. These data contain relatively large uncertainties in the glacial history and erosion on the Earth's surface. In this study, we use a highresolution present-day topography map, combined with the reconstructed bed deformation from a GIA model.
The reference PD topography is the RTopo-2 dataset by Schaffer et al. (2016), which is a compilation of various topography and bathymetry data sets on a 30 arcsec grid. Lake bathymetry, except the Great Lakes and Lake St. Clair, is missing from this data set. Consequently, other lakes cannot properly be reconstructed. For further processing, the data are remapped onto a lower resolution projected grid: a Lambert azimuthal equal area projection centred at 60 • N and 94 • W with a resolution of 5 km. By simply interpolating, sub-grid information about possible spillways would be lost. Therefore, a minimum-filter is applied to the high-resolved grid and the result is subsequently interpolated onto the projected grid. The minimum-filter assigns each cell the height of the lowest cell within a certain distance. It should be noted that this filtered map is only used to determine the lake basins, which are later transferred back to the unfiltered map.

Ice sheet and bed deformation
To represent the palaeogeography of the North American continent, the NAICE ice sheet reconstruction (Gowan et al., 2016b) is used. The ice sheets are determined using the ICESHEET model (Gowan et al., 2016a), which takes ice margin reconstructions and a model of basal shear stress as inputs and uses simplified ice physics to calculate ice surface elevation. A temporally and spatially varying basal shear stress model is tuned to create a glacial history that fits best with geological and geophysical measurements, such as GIA uplift rates, relative sea level and glacial lake strandline tilts. Gowan et al. (2016b) focused mainly on the western LIS. In this region the ice margins are inferred from the minimum timing of retreat given in Gowan (2013). For the eastern LIS and the Cordilleran Ice Sheet (CIS) the ice margins from Dyke (2004) are assumed. The reconstructions of the Innuitian and the Greenland Ice sheet were also obtained using the ICESHEET model by Khosravi (2017). The output of the NAICE model contains maps of the ice thickness and the bed deformation, which are combined with modern topography to produce the palaeogeography.
The proper reconstruction of post-LGM ice distribution and GIA requires estimates of the pre-LGM glacial history which is generally only poorly constrained (Gowan et al., 2016b). The NAICE reconstructions could, for example, be improved by using a more recent compilation of pre-LGM ice constraints (e.g. Batchelor et al., 2019).

LakeCC Algorithm
Depressions in a DEM act as a sink for many algorithms of flow routing. One possibility to overcome this issue is to remove all depressions from the DEM by filling them (Barnes et al., 2014;Barnes, 2016). Other, even more elaborate methods, create a drainage hierarchy (Barnes et al., 2019). The focus of this study is the identification of lake basins, without including advanced flow routing techniques. The LakeCC model is developed only for this purpose. For some other lake filling methods (e.g. Berends and van de Wal, 2016) a cell within each lake's basin has to be picked to determine its extent and minimum spillway. This is not required in our model. It assesses every lake basin of the domain, which is advantageous to get an overview of a whole map, without selecting individual lakes. The LakeCC model iterates over a set of increasing water levels, similar to the approach described in Wu et al. (2015) and Wu and Lane (2016).
The LakeCC algorithm requires a mask that contains information about sinks (i.e. ocean cells or domain margin) of the map. A simple way to do this would be to just set every cell below sea level to be an ocean cell. However, this would also label cells as ocean that are below sea level but inside the continent, and therefore prevent a lake from forming. A more advanced method, called SeaLevelCC, is also implemented into the LakeCC tool box.
Both methods, LakeCC and SeaLevelCC, are based on a simple 4-neighbour connected-component (CC) labelling algorithm, which is efficiently implemented using run-length encoding and union-find strategies (see Fig. 1). It was adapted from a tool written by Khroulev (2017) for identifying icebergs from a mask within the Parallel Ice Sheet Model (PISM) (the PISM authors, 2015). This underlying code iterates through the map and adds cells that satisfy the flooding criterion (Equation (1)) for a given water level h either to an adjacent patch or, if not present, to a new patch.
Here H is the bedrock elevation, T is the ice thickness and ρ ice and ρ water are the densities of ice and water, respectively. Additionally, if Equation (1)  the flooded cell is marked as a sink, the corresponding patch and all associated cells become a sink (see Fig. 1b, step II). The SeaLevelCC method requires a mask, with cells guaranteed to be ocean marked as sink. If this information is not available, the entire domain's margin can simply be set to sink before applying the CC method at sea level. By doing so the formation of inland ocean basins is effectively inhibited, as ocean patches need to be below sea level and connected to the domain's margin. After the land/sea mask is determined, the LakeCC algorithm is applied, which is shown in Fig. 2. It applies the CC method for a set of increasing water levels h. The resulting patches that are not labelled as sink represent a lake basin at the respective level. These patches grow in size as the level rises. Until a patch merges with another one that is a sink, i.e. the lake overflows, the current level h is stored for all associated cells (see blue box in Fig. 2). After this is done for various levels between h min and h max , the map contains an estimation of the basins' maximum fill levels within this range. The maximum uncertainty, apart from uncertainties originating from the discretised input topography, depends on the spacing dh between successive levels. For each lake, it is further checked if at least one cell of the lake is ice-free, to prevent sub-glacial lakes from forming.

Application
Fig. 2 depicts the application of the LakeCC algorithm as a flowchart diagram. The topography map is prepared as described above. Different filter diameters were tested and finally a value of 10 km was chosen as it makes the best compromise between draining under-resolved valleys and retaining realistic lake reconstructions. The domain used in this study is 7755 km x 5955 km, which, at 5 km resolution, is represented by a 1551×1191 rectangular grid. Successive water levels are evenly distributed between 0 and 1800 m by dh = 0.2 m.
The minimum-filtered topography map is assumed to retain information about possible spillways as needed for the filling algorithm. For the further steps, a smooth approximation of the actual bed is preferred. The lake levels obtained from the filtered map are transferred back onto the unfiltered map. This is possible since the filtered map represents an over-deepened version of the unfiltered map and therefore lake levels tend to be lower and thus still be valid on the unfiltered map.
The runtime of the LakeCC algorithm depends on the topography and ice configuration. It is shortest when large parts of the map are ice-covered and only a few lakes appear, and longest if many lakes exist in the domain. In our study, iteration over 9000 different levels took between 250 and 330 s on a 2.4 GHz Intel Xeon processor.

Post-processing
To get an estimate of the area, volume and maximum depth of certain lakes, another simple algorithm is applied to the output of the LakeCC model. It iterates over the grid and collects this data for all lakes.
In the next step, the map is divided into basins according to the slope of the topography. Each basin is assumed to overflow at its lowest boundary cell into a neighbouring basin. The drainage route is followed until a basin is reached that is marked as a sink (i.e. ocean or domain boundary). To prevent getting stuck in an infinite loop, it checks whether a basin has already been visited and marks the involved basins if necessary. In a subsequent iteration, all labelled basins are checked for the next highest spillway that escapes the loop. After iterating until no loops are left, each basin can be assigned a drainage route into the ocean.

Results
In this section, the lakes calculated for the different ice sheet reconstructions between the LGM and PD are presented. Some regions of interest are selected and highlighted in the following subsections. These regions are shown in the overview map (Fig. 3) as black rectangles. Overview maps of the ice-and lake reconstructions on the North American continent between 21 ka BP (Fig. S2A) and PD (Fig. S2V) are provided in S2. Although pre-LGM ice reconstructions are available from the NAICE model (Gowan et al., 2016b) these are not examined, since the ice margins are only poorly constrained. Table S1 gives the lake surface area, volume, water level, maximum depth and the sink (ocean basin where the lake drains to) for all identified lakes for all time slices.

Great Lakes
During the LGM the southern ice margin of the LIS was south of the Great Lakes, whose basins were fully ice-covered. As ice retreated, parts of the basins were exposed and became proglacial lakes. The reconstruction of glacial lakes in the Great Lakes' basins, including Lakes Maumee, Chicago, Algonquin, Duluth, Chippewa and Stanley, is in good agreement with the literature (Taylor, 1915;Larsen, 1988;Lewis and Anderson, 1989;Lewis et al., 1994;Breckenridge, 2013). Fig. 4 shows some snapshots of this region for different times. Fig. 5 shows that there is a good match of the PD lake reconstruction with the actual shorelines. Table 1 presents a more quantitative comparison of the lake reconstruction Major PD lakes are shown as thin black outlines with a black label: At -Athabasca; Er -Erie; GB -Great Bear Lake; GSa -Great Salt Lake; GSl -Great Slave Lake; Hu -Huron; LW -Lake of the Woods; Mb -Manitoba; Mi -Michigan; Ni -Nippigon; On -Ontario; Rd -Reindeer; Su -Superior; Wi -Winnipeg; Wp -Winnipegosis. The red outlines represent the basins of major palaeolakes discussed in the paper: Ag -Agassiz (Teller and Leverington, 2004); Bv -Bonneville (Chen and Maloof, 2017); Lh -Lahontan (Reheis, 1999b); MC -McConnell (Smith, 1994); MK -McKenzie (Smith, 1992); Ms -Missoula (Hanson et al., 2012); Pe -Peace (Mathews, 1980;Hickin et al., 2015). Areas where surficial sediments (Fullerton et al., 2003;Soller et al., 2009 (LakeCC) and measurements (CCGL). Since Lakes Huron and Michigan form a connected hydrological system, they are treated as a single lake. For all lakes, the reconstructed water level is close, but below the measured level, except for Lake Huron-Michigan, which is slightly overestimated. The same is true for the lake volume. Deviations in volume and surface area are between 1-4%, except for the relatively shallow Lake Erie, where the difference in volume is~80 km 3 (~16% different). The total area and volume of the Great Lakes are underestimated by 2.2% and 1.6%, respectively. Lake St. Clair is located north of Lake Erie's western end. Due to its size it is usually not counted among the Great Lakes and is not highlighted by its outline here. However, the reconstruction fits well with the actual lake geometry (see Fig. 5).

Great Basin
The Great Basin is located between the Rocky Mountains and the Sierra Nevada, mainly in Nevada, Utah and California (see rectangle 4 in Fig. 3). Although the ice margin never extended this far south during the time span of interest, this region is treated here because the mountainous terrain has several lake basins.
The lake reconstructions show the formation of huge lakes in this area for every time slice; see Fig. 6a and overview maps ( Fig. S2A-V). In Fig. 6a the red outlines mark the maximum extent of Lake Bonneville (Chen and Maloof, 2017), the predecessor of the Great Salt Lake, and Lake Lahontan (Reheis, 1999b), which covered the Black Rock Desert. These outlines match quite well with the reconstructed lake extents. Apart from these two lakes, the reconstruction also shows other filled valleys: the long lake parallel to the border between Nevada and California filling up Death Valley coincides with Lake Manly (compare with Blackwelder (1933) and the lacustrine sediments shown in Fig. 3). Also, the canyons cut by the Colorado River in southeastern Utah are filled up.

Central North
At the LGM the southern margins of the Laurentide and Cordilleran ice sheets were located in central Montana (Dyke, 2004), along which several small lakes formed, as can be seen in Fig. 6b (this region corresponds to rectangle 2 in Fig. 3). The red outline shows the maximum extent of Lake Missoula (Pardee, 1910;Hanson et al., 2012). This lake, however, is not captured in the reconstruction; only small separate lakes fill up in this basin. The locations of the small lakes east of Lake Missoula along the ice margin are in good    Fig. 7a shows an overview of the central north of North America at 15 ka BP. In this time slice, the LIS and CIS are about to separate. The reconstruction predicts the formation of huge lakes with a total volume of about 118 000 km 3 of freshwater in areas that just deglaciated in this time slice. For these lakes, however, no geological evidence is available. Due to lack of a name, these are called X and Y here. On the southern end of the LIS, the ice retreated far enough to expose parts of Lake Agassiz's basin. This newly formed lake fits quite well into the proposed margin of Teller and Leverington (2004). The small lake west of Lake Agassiz coincides well with Lake Souris (Kehew and Teller, 1994;Sun and Teller, 1997).
In Fig. 7b, the reconstruction of this region for 13 ka BP is shown. In the northwest, the retreating LIS exposed the Great Bear Lake basin, which contains glacial Lake McConnell. Where lake Y was previously, two much smaller lakes are found: Lake Meadow and Lake Saskatchewan (Christiansen, 1979;Kehew and Teller, 1994). In northwestern Alberta, Lake Peace (Hickin et al., 2015) formed. Table 2 contains surface area, volume, maximum depth, and the sink of the major lakes found in this region for the time between 15 ka BP and 9 ka BP. Lake Agassiz is found to drain towards the south until 12 ka BP (see Fig. 7c). In this time slice, the retreating ice margin exposed the basin of Lake McKenzie and parts of Lake McConnell's basin. Lake McConnell fits nicely into the basin reported by Smith (1994), but Lake McKenzie does not fill up except for a few separate cells. Fig. 7d shows the reconstruction for 11 ka BP. The LIS retreated so far that a new drainage route for Lake Agassiz opened towards the Arctic ( Table 2). The lake followed the ice margin and thus nearly doubled in size. It stays inside its proposed margin of maximum extent (Teller and Leverington, 2004), as does Lake McConnell (Smith, 1994).
At 10 ka BP (Fig. 7e) the ice margin in the west has not changed much compared with the previous time slice, but the margin in the south and southwest has retreated several hundreds of kilometres and the overall ice volume decreased. Lake Agassiz reaches its maximum size for the snapshots analysed here by almost doubling its volume once again (see Table 2). In this configuration, it is draining towards the south. Lake McConnell's size decreased, and is split into separate basins.
In Fig. 7f (9 ka BP), the southern margin of the LIS rapidly retreated northward. This opens a new basin in the east for Lake Ojibway, which fits into the maximum outline of the Lake Agassiz/Ojibway basin. It appears twice in Table 2 because the two parts in Fig. 7f are not connected. Lake Agassiz is dramatically reduced in extent. Both lakes drain via the St. Lawrence River into the North Atlantic. This is the last time slice with major proglacial lakes appearing. At 8 ka BP, the LIS over Hudson Bay collapsed and overall dramatically decreased in size.

Testing the LakeCC method
To validate the LakeCC model, the PD reconstructions of the Great Lakes are discussed in more detail. The reconstructed lake extent fits nicely with the actual shorelines (see Fig. 5) and also compares well with measurements (Coordinating Committee on Great Lakes Basic Hydraulic and Hydrologic Data, 1977) shown in Table 1. It should be noted that part of the difference is because different lake bathymetry data sets were used. However, we assume that this contribution is minor. The approximately 2% deviation in total area and volume are assumed to be due to the relatively coarse 5 km resolution. Deviation in water level is related to resolution that prevents the drainage route from being properly resolved. This issue results in small shallow lakes forming and can only be overcome by applying a minimum-filter to the high-resolution map, as was described above. This yields good results: major lakes are retained, while most unresolved valleys are drained. However, this comes at the cost of also draining lakes in narrow valleys that were actually existing, such as Lake Missoula (Fig. 6b) or Lake McKenzie (Fig 7b-f). Fig. 6a and the overview maps ( Fig. S2A-V) show that the lakes in southeastern Utah experienced only minor changes in time due to GIA. Topography reconstructions in this area are therefore similar to the PD reference topography. As the canyons are drained by the Colorado River, we assume that the lakes in our reconstruction are artefacts originating from the coarse resolution.
The LakeCC tool basically only determines lake basins. How these are interpreted is up to the user. In this study, these lake basins are assumed to be filled to the top, for simplicity, but also because of a lack of any further information, such as climate and a temporally higher resolved glacial history. This simple assumption seems to be valid along the ice margins at times of glacial retreat. Examples for this issue being a problem in this study are the presence of Lake Bonneville (the predecessor of the Great Salt Lake), Lake Lahontan (in the Black Rock Desert) and Lake Manly (in Death Valley) at each time slice (see above). These are located in the Great Basin, which at present is the largest endorheic basin in North America and is characterised by an arid climate. The climatic conditions in this area were different in the past, as strandlines and lake sediments provide evidence of large ancient lakes (Chen and Maloof, 2017). These are assumed to have evaporated as the climate became drier (Reheis, 1999a). The reconstructed lakes, however, fit well within the outlines of maximum extent (see Fig. 6a).
Despite these uncertainties, the reconstructions are quite valuable in the context of this study, as we were more interested in the formation of the major lakes along the ice margin. A more elaborate study design would require more detailed information that is not available for such palaeo-applications.

Validation of the ice model
For this study, ice sheet and GIA reconstructions produced with the ICESHEET model (Gowan et al., 2016a,b) were used to reconstruct lakes forming along the ice margin. This ice model was tuned to match various observations, such as GIA uplift rates, proglacial lake strandline tilt and relative sea level. One application of the present study is as a further indicator for the quality of the glacial reconstruction. As lake formation is quite sensitive to the ice margin location, the LakeCC model could be used to fine-tune ice margins in places where lakes form for which no geological evidence exists. Also, the water load should be included in the calculation of the GIA. This circular application is not done within this study; the lakes are rather used as a kind of benchmark to assess the glacial history of the ice sheet reconstruction. One of those lakes without geological evidence is repeatedly forming between the LIS and CIS at around 60 • N (see Fig.  7a), which, due to lack of a name, is called X in Table S1. It is present at the LGM (21 ka BP) but disappears again as the ice sheets merge at the beginning of the deglaciation. At 16 ka BP, the ice covering its basin has become thin enough that it fills again. The reason is that the Mackenzie River, which drains this region today, is fully glaciated. The lake disappears as this drainage route is open after 14 ka BP. Adjusting the ice margin of the LIS in this region a few kilometres to the east might solve this issue for 14 ka BP. For the earlier time slices this shift of the margin would be quite large. The reason why this basin did not fill up might simply be that the climatic conditions in this region were too cold.
The other lake is called Y in Table S1, and only forms between 15 and 14 ka BP at the southwestern margin of the LIS.
It is assumed that by shifting the ice margin southeast of this lake slightly to the northeast could drop the lake level about 200 m, which would almost drain the entire lake, leaving only small separate lakes, such as an early stage of Lake Peace at 14 ka BP in central Alberta (compare with Hickin et al. (2015) and the surficial sediments shown in Fig. 3). Additionally, the model was applied to the ICE-6G_C reconstructions (Peltier et al., 2015). As results are sensitive to the ice margin, the roughly resolved ice sheet was interpolated onto a higher resolved grid and cropped according to the ice margins used for the creation of ICE-6G_C (Dyke et al., 2002(Dyke et al., , 2003. Compared with the reconstructions of the NAICE model (Gowan et al., 2016b) ICE-6G_C has more ice. This, in addition to the fact that a different Earth model structure was used, explains the deeper depression of the surface, which can be seen in the lake reconstructions. Fig. 8 shows the lake reconstruction for 12 ka BP; all other time slices can be found in S2 (Fig. S2A-V). The southern tip of Lake Agassiz in Fig. 8 is north of the Moorhead Low margin, that occurred at about that time (Breckenridge, 2015). The western margin merges with Lake McConnell and Lake McKenzie to form an enormous lake. The basins of the Great Lakes appear strongly tilted northward, with only small lakes forming in the basins of Lake Michigan and Lake Erie. Lake Ontario merged with the Champlain Sea.

Drainage routes
Input of freshwater into the ocean at certain key regions can disrupt the thermohaline circulation and thus impact the global climate system (e.g. Manabe and Stouffer, 1997;  Lohmann and Schulz, 2000;Kageyama et al., 2009;Condron and Winsor, 2012). Broecker et al. (1989) proposed that the rerouting of LIS meltwater away from the southern spillway via the Mississippi River towards the east via the St. Lawrence Stream might have caused the onset of the Younger Dryas (YD). The diversion of drainage is consistent with a δ 18 O anomaly measured in sediment cores from the Gulf of Mexico, which was caused by a sudden lack of freshwater input at the beginning of the YD (~12.9 ka BP) (Wickert et al., 2013). However, there is still debate as to whether this rerouting of Lake Agassiz overflow happened towards the east (Carlson et al., 2007;Wickert et al., 2013;Leydet et al., 2018) or towards the Arctic via the Mackenzie River (Murton et al., 2010;Keigwin et al., 2018). These different drainage routes are marked as green arrows in Fig. 3. Applying the post-processing described above to our LakeCC reconstructions we estimated the drainage routes and affected ocean basins for Lake Agassiz (see Table 2). From the first appearance of Lake Agassiz at 15 ka BP it overflows at its southern end and drains via the Mississippi River into the   Figure 8. This map gives an overview of the lake reconstructions in North America at 12 ka BP using the ICE-6G_C ice and palaeogeographic reconstruction (Peltier et al., 2015) [Color figure can be viewed at wileyonlinelibrary.com].
Gulf of Mexico. At the beginning of the YD at 13 ka BP (Fig. 7b) the northern half of Saskatchewan is still covered by thick ice, as is the northern half of the Great Lakes (see Fig. S2Ia). This prevents the lake draining to the north and east. The southern drainage route in our reconstructions is active until 12 ka BP, after the onset of the YD. These results neither support, nor refute any of these drainage scenarios, as the lake outlets depend strongly on an accurate reconstruction of the ice margin. The southeastern LIS margin used for the NAICE reconstruction (Gowan et al., 2016b) resembles the ice margin by Dyke (2004), i.e. recent discoveries about the timing of ice retreats (e.g. Leydet et al., 2018) are not accounted for. Furthermore, due to the coarse temporal resolution of the ice sheet reconstruction (1 ka), advances and retreats of ice lobes are not resolved. The ice margin blocking the northwestern drainage route does not necessarily contradict the freshening of the Beaufort Sea at the beginning of the YD, as reported by Keigwin et al. (2018). The source of this signal can be attributed to the rapid deglaciation of the northwestern LIS Peltier, 2005, 2006) and drainage events of other major proglacial lakes in this region, such as Lake McConnell or Lake Peace (Wickert, 2014(Wickert, , 2016. Breckenridge (2015) dated the palaeostrandline of Lake Agassiz that is associated with the opening of a lower spillway towards the north to 12.18 ± 0.48 ka BP. This happens in our reconstructions at 11 ka BP, when the northwestern part of the LIS has retreated far enough to open a new gateway to the Mackenzie River (see Fig. 7d). The outlets to the north and south do not differ much in height, as the southwestern shoreline of Lake Agassiz barely changed compared with the previous time slice (Fig. 7c). At 10 ka BP (Fig. 7e) the northward drainage route has closed again, possibly due to GIA, and drainage is again via the Mississippi River. The re-activation of the Mississippi drainage route is consistent with Wickert (2014). As already mentioned before, the exact timing is sensitive to the ice margin positions, which are often poorly constrained. Especially dynamic processes, such as advances and retreats of ice lobes, are only poorly covered.
Before Lake Agassiz finally drains around 8 ka BP as the LIS over Hudson Bay breaks apart, a new drainage route via Lake Ojibway to the St. Lawrence River appears (9 ka BP, Fig. 7f). In our reconstruction, the opening of this spillway yields an extreme drop in lake level. The diversion of flow away from Lake Superior towards Lake Ojibway can be found in sedimentary records and is dated to ∼9.04 ka BP (Breckenridge et al., 2004(Breckenridge et al., , 2012Wickert, 2014).

Conclusions
In this study we have presented a tool to efficiently locate closed basins for a given palaeotopography configuration. After pre-processing the topography data set, we are able to retrieve lake reconstructions that, despite the relatively coarse grid, are in good agreement with geological data. For PD, the reconstructed Great Lakes lake volume, surface area and water levels agree well with observations. The LakeCC model was also applied to palaeoreconstructions of ice and GIA from the NAICE (Gowan et al., 2016b) and ICE-6G (Peltier et al., 2015) models. The results show the evolution of major glacial lakes, such as Lake Agassiz, Lake McConnell and the various predecessors of the Great Lakes. The discrepancy between reconstruction and geological data was used to infer places where the ice margin and volume might need adjustments.
The LakeCC model can be used to improve the ice sheet and topography reconstructions by applying adjustments of the ice margin and ice volume, and by adding the water load of the lakes into the calculation of the GIA. Furthermore, the model results give an idea about deglacial meltwater distribution, which is an important factor in the global climate system. In a future study we will couple the LakeCC model to an ice sheet model and investigate the impact of lakes on ice dynamics.

Supporting information
Additional supporting information may be found in the online version of this article at the publisher's web-site. S1 Table containing area, volume, lake level, maximum depth,and sink for all identified lakes S2 Overview maps of the lake reconstructions (for both, NAICE and ICE-6G) in North America for different time slices