Reconstruction of North American drainage basins and river discharge since the Last Glacial Maximum

Over the last glacial cycle, ice sheets and the resultant glacial isostatic adjustment (GIA) rearranged river systems. As these riverine threads that tied the ice sheets to the sea were stretched, severed, and restructured, they also shrank and swelled with the pulse of meltwater inputs and time-varying drainage basin areas, and sometimes delivered enough meltwater to the oceans in the right places to influence global climate. Here I present a general method to compute past river flow paths, drainage basin geometries, and river discharges, by combining models of past ice sheets, glacial isostatic adjustment, and climate. The result is a time series of synthetic paleohydrographs and drainage basin maps from the Last Glacial Maximum to present for nine major drainage basins – the Mississippi, Rio Grande, Colorado, Columbia, Mackenzie, Hudson Bay, Saint Lawrence, Hudson, and Susquehanna/Chesapeake Bay. These are based on five published reconstructions of the North American ice sheets. I compare these maps with drainage reconstructions and discharge histories based on a review of observational evidence, including river deposits and terraces, isotopic records, mineral provenance markers, glacial moraine histories, and evidence of ice stream and tunnel valley flow directions. The sharp boundaries of the reconstructed past drainage basins complement the flexurally smoothed GIA signal that is more often used to validate ice-sheet reconstructions, and provide a complementary framework to reduce nonuniqueness in model reconstructions of the North American ice-sheet complex.


15
At the time of Last Glacial Maximum ice extent, ca. 30 ka to 19.5 ka Lambeck et al., 2014), North American drainage configurations and river discharge were dramatically different from present. Ice advance rerouted rivers (e.g., Ridge, 1997;Curry, 1998;Blumentritt et al., 2009), and ice streams fed temporarily-enlarged drainage basins (Dyke and Prest, 1987;Patterson, 1998;Margold et al., 2014Margold et al., , 2015. The growth of the ice-sheets changed atmospheric circulation by 20 generating up to ∼4 km of high-albedo ice-surface topography (Kutzbach and Wright, 1985;Ullman et al., 2014) that influenced global climate and therefore patterns of drainage-basin-scale water balance and river discharge to the oceans.
Reconstructions of past ice-sheet thickness have proliferated (e.g., Tushingham and Peltier, 1991;Lambeck et al., 2002;Peltier, 2004;Tarasov and Peltier, 2006;Tarasov et al., 2012;Gregoire et al., 25 2012; Argus et al., 2014;Peltier et al., 2015), but are evaluated using limited (Tarasov and Peltier,1 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-8, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 February 2016 c Author(s) 2016. CC-BY 3.0 License. 2006) to no (all others) geologic evidence for past drainage patterns. Instead, they are tested against terminal moraine positions and glacial isostatic adjustment, with the latter being a response that is spread across hundreds of to ∼1000 kilometers due to the high flexural rigidity of the lithosphere, especially in the Canadian Shield provence that was covered by the ice-sheet (Kirby and Swain, 2009;30 Tesauro et al., 2012). Flow-routing calculations (e.g., Metz et al., 2011;Schwanghart and Scherler, 2014) are deterministic and can be used to produce time-series of drainage pathways predicted by each reconstructed ice-sheet and its associated glacial isostatic adjustment (GIA) history, forming the necessary links between the reconstructed ice-sheet and measurable geomorphic and sedimentary records. Furthermore, the discrete boundaries of drainage basins are sensitive to ice geometry in 35 a way that is distinct from the broad GIA response, and therefore can provide a partly-independent check on the accuracy of any given ice-sheet reconstruction (Wickert et al., 2013).
On a global scale, climate change during deglaciation may have been profoundly affected by meltwater delivery to the oceans, especially at the time of the Younger Dryas (e.g., Broecker et al., 1989), but is highly sensitive to the locations of meltwater inputs to the oceans (e.g., Tarasov and Peltier,40 2005; Murton et al., 2010;Condron and Winsor, 2012). Geologic data provide a timeline of meltwater delivery to one basin or another (e.g., Ridge, 1997;Knox, 2007;Rittenour et al., 2007;Murton et al., 2010;Williams et al., 2012;Breckenridge, 2015), but do not provide estimates of volume.
Furthermore, the relationships between geological data in different parts of the continent must be connected by conservation of ice-sheet mass and glaciological continuity of ice flow. Continental-45 scale meltwater routing calculations that connect data and models address both of these concerns by using quantified ice-sheet thickness and ensuring that water that is not sent to one basin is sent to another, in essence serving as a "connective tissue" to ensure quantitative self-consistency.
In general, there exists a lack of recognition of the importance of Pleistocene drainage rearrangement on river systems. In spite of the broad recognition that river systems and drainage basins have 50 changed significantly since the LGM (Dyke and Prest, 1987;Wright, 1987;Teller, 1990a, b;Marshall and Clarke, 1999;Patterson, 1997Patterson, , 1998Licciardi et al., 1999;Overeem et al., 2005;Tarasov and Peltier, 2006;Wickert, 2014;Margold et al., 2014Margold et al., , 2015, many studies do not include any welldefined picture of drainage basin evolution. This results in representations of the modern drainage network wholly or partially in place of glacial-stage drainage basins (e.g., Sionneau et al., 2008;55 Montero- Serrano et al., 2009;Sionneau et al., 2010;Lewis and Teller, 2006;Tripsanas et al., 2007) that help to propagate a lack of consciousness about the continental-scale hydrologic changes that occurred, even in cases when the study acknowledges in the text that past drainage pathways were different. Model-based studies likewise either test a range of possible drainage basin areas (Overeem et al., 2005) or assume that the rivers have remained unchanged (Roberts et al., 2012). In the former 60 case, this can lead to less-precise calculations of sediment transport and deposition. In the latter case, this violates the underlying assumptions used to compute erosion and infer spatially-distributed uplift. All of these scenarios highlight the need for an improved community-wide understanding of the The drainage basin and paleohydrologic reconstructions are a combination of four inputs: (1) a digital elevation model (DEM) with high enough resolution to resolve the valleys of major river systems, (2) a prescribed ice-sheet history, (3) a model of the GIA response to that ice history, and 100 (4) a past water balance based on general circulation model (GCM) outputs that have been offset to fit modern data at the present time step (Fig. 1. These factors together allow us to generate synthetic river discharge histories that can be compared with geologic data. This analysis produced three major end products. The first is a spatially-distributed map of estimated mean river discharges since the Last Glacial Maximum. The second is a gridded product of 105 river discharges to the sea, designed for coupling with general circulation models to investigate the feedbacks between ice melt, ocean circulation, and climate (see Ivanović et al., 2014). The third is a set of major drainage basin extents and a time-series of their discharges into the ocean.
The drainage analyses were performed using GRASS GIS (Neteler et al., 2012;GRASS Development Team, 2012), a free and open-source GIS software that has efficient data management, 110 computational, and hydrologic tools. GRASS GIS is fully scriptable using a number of languages.
The Python source code for these analyses has been released under the GNU General Public License (GPL) version 3 and is made available to the public at the University of Minnesota Earth-surface GitHub organizational repository, at https://github.com/umn-earth-surface/ice-age-rivers.  Water and ice flow downslope following a routing surface, Z r , that is given by: Z b is the reference (i.e., modern) land-surface elevation (subscript "b" stands for "modern glacier bed" or "modern bare (i.e., ice-free) ground"), H i is the ice thickness, and ∆Z b represents changes between the past and modern land-surface (glacier-bed) elevation. Together, Z b + ∆Z b give the 120 contemporaneous glacier bed elevation field. When this (Z m,b + ∆Z b ) surface is summed with ice thickness, H i , the resulting routing surface, Z r , combines ice-free land-surface elevations and the elevations of the surface of the ice-sheet. This ice-sheet surface is chosen to drive the flow-routing calculations because this is what drives ice flow (cf. Cuffey and Paterson, 2010), and therefore sets ice divides and meltwater pathways. If one wanted to compute subglacial water routing instead, an 125 alternate approach is to drive flow by the subglacial pressure gradient to simulate meltwater routing through subglacial streams and tunnel channels (e.g., Wright, 1973), and this could be accomplished by multiplying H i by ρ i /ρ w , where the latter are the densities of ice and water, respectively (e.g., Livingstone et al., 2015).  (van Hise and Leith, 1911;Dean and Phillips, 2011), indicate 145 that still-higher resolution elevation data will continue to improve flow-routing calculations near headwaters streams, a scale-dependence given by the relationship between channel geometry, river discharge, and drainage area (Leopold and Maddock, 1953).
Raster grid cell area, A c , is a function of latitude, and this is an important conversion factor for changing scalar quantities in cells, such as ice thickness, into a measure of volume. Using the 150 simplifying assumption that Earth is perfectly spherical, Here, R E is the mean radius of Earth, θ is again colatitude, and ϕ is again east-longitude (i.e., degrees east of Greenwich). δθ is cell north-south extent, and δϕ is cell east-west extent; in the simulations presented here, both of these extents equal 30 arcseconds. These areas are computed at a latitude λ 155 corresponding to the north-south midpoint of the cell, with the implicit assumption that the cells are small enough that the north-south difference in cell size can be approximated to be linear such that a straight mean accurately represents the whole. For the GEBCO_08 30-arcsecond grid that is used here, this approximation generates errors of <0.3%.

160
H i , ice-sheet thickness, is provided by pre-generated ice-sheet reconstructions. These are produced at a coarser resolution than the 30-arcsecond flow-routing grid, and therefore are interpolated using an iterative nearest-neighbor approach to remove stepwise discontinuities that would otherwise introduce artifacts in the flow-routing calculations.
I use five ice models, four of which were created based on geophysical data, and one of which 165 was generated using a numerical model of ice-sheet physics. In the former, ice-sheet thickness evolution was determined by inverting for global sea level data by combining (1) the total water mass contained within the ice-sheets, causing global mean sea level fall, and (2) calculations of the GIA process, which created deviations from the mean sea level history. Most notable among these deviations is that due to postglacial rebound, which is a key indicator of the presence of major past 170 ice-sheets. However, postglacial rebound rates are a function of both time-evolving ice-sheet thickness and solid Earth (lithosphere and mantle) rheology. That this one rate is a function of two only partly-constrained variables introduces an element of uncertainty into GIA-based inversions (Mitrovica, 1996). Ice physics, on the other hand, are sensitive to the thermal evolution of the ice-sheet, related thermomechanical effects on ice rheology, and the ice basal conditions-a large number of 175 unknowns (Tarasov and Peltier, 2006;Cuffey and Paterson, 2010). Therefore, I investigate both types of models.
ICE-3G (Tushingham and Peltier, 1991) was one of the earliest widely-used ice-sheet reconstructions. It was constructed to constrain ice thickness as much as possible with measurements of glacial isostatic adjustment from radiocarbon-dated sea level markers. While it fails to reproduce the full 180 magnitude of the observed LGM global sea level fall (Austermann et al., 2013;Lambeck et al., 2014), it provides a counterpoint to more recent models that include a much more massive Laurentide ice-sheet (cf. Lambeck et al., 2002;Peltier, 2004).
ICE-5G (Peltier, 2004) has been often considered the "industry standard" global ice-sheet and GIA reconstruction. It was also built by combining geological data and geophysical inversions, and 185 is an update to ICE-3G. It contains a very thick north-south-oriented ice dome over western Canada (Argus and Peltier, 2010) (Wickert et al., 2013), and in any case, seems unlikely as the ice-sheet would slide and thin over the soft Western Interior Seaway sediments east of the Canadian 190 Rockies. Nevertheless, ICE-5G remains widely-used and matches many observations (Peltier, 2004).
ICE-6G (here, shorthand for ICE-6G_C) (Argus et al., 2014;Peltier et al., 2015) is the successor to ICE-5G. It has been shown to improve agreement with present-day GIA measurements, and matches geological evidence for the locations of the ice domes and ice-sheet outlines (Patterson, 1998;Dyke et al., 2003;Dyke, 2004). As ICE-6G is so new, it has not been extensively tested outside of the 195 group that developed it (Argus and Peltier, 2010;Argus et al., 2014;Peltier et al., 2015). Lambeck et al. (2002) developed the Australian National University (ANU) ice-sheet model, also based on field data constraining sea level and GIA. This model differs from ICE-3G and ICE-5G in that it was developed regionally, and these regions were later combined into a more global picture of ice-sheet evolution.

200
To produce the G12 ice model, Gregoire et al. (2012) simulated North American ice-sheet evolution from the LGM to present using the Glimmer community ice-sheet model (Glimmer-CISM) (Rutt et al., 2009). Glimmer-CISM is a thermomechanical model of ice deformation and dynamics that uses the shallow ice approximation for its driving stresses. While the shallow ice approximation neglects higher-order (i.e., less-important) stresses, it generally closely approximates those models 205 that do include higher-order terms, and its computational efficiency makes it a viable option for long-term simulations (Rutt et al., 2009), such as the Gregoire et al. (2012) LGM-present simulation. One of its major features is a collapse of the ice saddle between the Laurentide and Cordilleran Ice Sheets at ca. 11.6 ka. Geological data, on the other hand, place the major phase of saddle collapse during the Bølling-Allerød warm period, ca. 14.5-12.9 ka Williams et al. (2012); Dyke et al.

210
(2003); Dyke (2004), and possibly beginning as early as ∼17.0 ka (John et al., 1996;Clague and James, 2002;Carlson and Clark, 2012). While Gregoire et al. (2012) shifted their chronology to compare their Laurentide-Cordilleran ice-saddle collapse to the geologic record of Meltwater Pulse 1A (Fairbanks, 1989;Deschamps et al., 2012), I use their model ages as direct geologic age. Most importantly, the ice dynamics contained within Glimmer-CISM provide it a set of tools to match 215 reality that place it in a different category than the GIA-based models.
2.1.3 Change in land-surface elevation: glacial isostatic adjustment ∆Z b , change in land-surface (glacier bed) elevation, is here provided by the combination of changes in global mean sea level (GMSL) and GIA. At the Last Glacial Maximum, GMSL was ∼135 m lower than present, but local relative sea level was highly variable due to the effects of GIA (Lambeck et al.,220 2014). GIA comprises deformational, gravitational, and rotational effects of ice-age loading. The deformational component is the local flexural isostatic response to changes in ice and water loading.
The gravitational component is in response to redistribution of mass (and hence, gravitational potential) across the surface of Earth, and this is amplified by the water self-gravitation feedback (e.g.,
Published: 17 February 2016 c Author(s) 2016. CC-BY 3.0 License. Mitrovica and Milne, 2003). The rotational component is also the response to changing near-surface 225 mass distributions that perturb the moment of inertia and produce true polar wander that causes the equatorial bulge to migrate, first in the form of seawater and later in the form of solid Earth material.
These three components feed back into one another, and therefore are solved numerically through iteration.
This gravitationally self-consistent sea level theory and its numerical implementation are de-230 scribed by Mitrovica and Milne (2003) and Kendall et al. (2005), respectively, who calculate the GIA effects of changing ice and water masses. The calculations are performed using a pseudospectral algorithm that, for this work, is truncated at spherical harmonic degree and order 256. This satisfies the Nyquist flexural wavenumber for the solid Earth models employed, thereby allowing the solutions to be interpolated to an arbitrarily high resolution.

235
ICE-3G, ICE-5G, ICE-6G, and ANU were each calibrated to match sea level data based on specific 1D models of solid Earth viscoelasticity (Tushingham and Peltier, 1991;Lambeck et al., 2002;Peltier, 2004;Argus et al., 2014;Peltier et al., 2015). I use their respective spherically symmetric models of Earth as a Maxwell body to simulate its response to these loads-VM1 for ICE-3G, VM2 for ICE-5G, VM5a for ICE-6G, and an unnamed model for ANU. G12 was built using a simplified 240 model of isostatic response to glacial loading, and was run only for North America. I pair it with VM2, and a posteriori enforce a total global ice-sheet volume (and therefore GMSL history) equal to that of ICE-5G. Among these solid Earth models, VM5a is largely a simplification of VM2, while the ANU solid Earth model has a higher viscosity contrast between the upper and lower mantle than VM2, placing it in better agreement with geophysical inferences of the upper mantle-lower mantle 245 viscosity contrast based on observations of the long-wavelength geoid (Hager and Richards, 1989).
I initiate the GIA modeling of G12 in equilibrium at the LGM, and ANU slightly before the LGM. ICE-5G includes hypothesized ice-sheet growth since 120 ka. This start point is during the last interglacial (MIS 5e), when GMSL peaked 5.5-7.5 meters above its present level (Dutton and Lambeck, 2012;Hay et al., 2014), though with some significant fluctuations (Hearty et al., 2007;250 Rohling et al., 2007;Dutton and Lambeck, 2012;Hay et al., 2014). ICE-3G was initiated at the last interglacial as well, with ice slowly grown in a simple pattern towards the LGM before starting to shrink.
Our flow-routing calculations neglect geomorphic change for three major reasons. First, significant changes Lake Agassiz outflow directions occurred as a precursor to spillway incision (Teller 255 et al., 2002), meaning that ice-sheet geometry and GIA are the primary drivers of drainage change.
Second, while process-based landscape evolution models exist (cf. Tucker and Whipple, 2002;Willgoose, 2005), the thresholded and highly nonlinear relationships between flow, erosion, sediment transport, and deposition, means that these models may predict a very wide range of potential landscape responses with sufficient uncertainty as to reduce their predictive capability. This would leave 260 changing the DEM's by hand as an option (as was done by Tarasov and Peltier, 2005), but consid- Each cell in the flow-routing grid domain can act as a source of water to continental drainage systems.
Water balance over each cell is a sum of precipitation gains, P , evapotranspiration losses, ET , and changes in water-equivalent thickness of the ice-sheet reconstruction, (dH i /dt) * (ρ i /ρ w ). When multiplied by the cell area, A c (Eq. 2), the result is an expression for volumetric incoming water

Global precipitation and evapotranspiration
Precipitation and evaporation inputs changed through time (Fig. 2). These changes were calculated by Liu et al. (2009) (Collins et al., 2006). This simulation, called TraCE-21K, was run at 3.5 degree resolution. I spherically interpolated (following Dierckx, 1993) the results to 0.25 degree (15 arcminute) resolution to match the assembled modern data products (below).
The low resolution CCSM3 contains systematic biases that produce a ∼1 mm/day positive precipitation anomaly in western North America and a ∼1 mm/day negative precipitation anomaly in 280 eastern North America (Yeager et al., 2006). To correct for this and other possible model biases, I assembled modern precipitation and evapotranspiration fields, computed the difference between each time-step of TraCE-21K and its modern time-step, and summed the modern compilations with the computed TraCE-21K difference from modern. I could find no global evapotranspiration product, and global precipitation products were coarse, often at 2.5-degree resolution (∼280-km cells at the 285 equator) (Xie and Arkin, 1997;Chen et al., 2002;Adler et al., 2003). Therefore, I generated my own global data product compilations as mean rates from 2000 to 2004, for consistency with prior work (Wickert et al., 2013), at 0.25 • resolution.
I derived the precipitation data by stitching together three data products. The CRU TS3.21 0.5 • precipitation product provided precipitation over the global continents except Antarctica (expanded 290 from the work of Harris et al., 2014), and this was smoothed and spherically interpolated (following Dierckx, 1993) to 0.25 • resolution. The HOAPS-3 satellite-derived precipitation data product (0.25 • resolution) was used for the oceans north of the Antarctic region (Andersson et al., 2010). Gaps of 2 • or less between these data sets were filled by an interpolation of them. The remainder of the globe, Our modern evapotranspiration rates for the continents (except Antarctica) are based on calculations using the MODerate-resolution Imaging Spectroradiometer (MODIS) instruments that are mounted on the Aqua and Terra satellites (Mu et al., 2007(Mu et al., , 2011, available at 30-arcsecond resolution. The oceanic evaporation data are from HOAPS-3 and at 0.25 • resolution (Andersson et al., present. These differences were calculated with the TraCE-21K continuous LGM to present paleoclimate run (Liu et al., 2009;He, 2011) of the CCSM3 GCM (Collins et al., 2006). They are combined with ice mass balance to produce water inputs for computed paleodischarges

Ice melt and central differencing
Ice-sheet inputs to runoff are determined by differencing ice-sheet thicknesses at adjacent time steps for each of the ice models. These differences in ice-sheet height are converted to water discharges by (1) multiplying them by the DEM cell size (∝ cos(latitude)), (2) multiplying them by 0.917 as a 310 conversion factor between ice and water density, and (3) dividing them by the time elapsed (Eq. 3).
Flow routing and drainage basins can be calculated only on time steps when the ice-sheet and GIA models provide surface topography (Section 2.3, below), but the ice-sheet contribution to runoff requires ice-sheet thickness to be differenced, and thus should be most valid for the times halfway between the ice thickness time-steps. This requires a choice of whether to allow numerical diffu-315 sion in space (i.e., flow-routing) or time (i.e., melt rate). Melt rates change more gradually than near-instantaneous drainage rerouting events, and these rerouting events are thought to be critical to understanding climate impacts of ice melt (e.g., Broecker et al., 1989;Tarasov and Peltier, 2006;Condron and Winsor, 2012). Therefore, I prefer numerical diffusion in time, and allow the gridded melt rate at each ice-sheet reconstruction time-step to equal the averages of the melt rates calcu-320 lated between this time step and those before and after it. This is the same as a central-difference calculation when time-steps are uniform.

Neglected terms: groundwater and lakes
While groundwater is a significant reservoir that affects river discharge, I assume that over the millennial time scales of interest to us, groundwater discharge will be in equilibrium with atmo-325 spheric inputs. Although ice-sheets can strongly influence groundwater discharge, modeling work by Lemieux et al. (2008) shows that the maximum net deglacial groundwater discharge across the whole North American continent is only ∼2.5% the modern Mississippi River discharge, and that this spike was short-lived.
I do not include proglacial lakes in the continental water storage because the ice-sheet reconstruc-330 tions discussed here (Tushingham and Peltier, 1991;Lambeck et al., 2002;Peltier, 2004;Gregoire et al., 2012;Peltier et al., 2015) do not. This means that these models implicitly include the lakewater volume within the ice mass. This deficiency affects the ability of the models to match the observed ice-sheet geometry, as ice must be used in the models to represent the surface loads of these large lakes in the geologic past, and removes the potential effects of the lakes in driving ice-sheet retreat in mechanistic models. In this work, water is simply routed across proglacial-lake-producing depressions (see Section 2.3).

Continental flow routing and accumulation
After defining the amount of flow each cell contributes to runoff, I used the time-evolving grids of surface elevation to compute drainage patterns across North America. I first specified that the 340 drainage outlets must occur along the coast at that time-step. I then routed flow down the Z r (Eq. 1) surface while computing accumulated discharge as (1)

Defining continent, ocean, and shore
Shorelines and topographic elevations changed dramatically since the Last Glacial Maximum, thus requiring redefinition of the coastline on which river mouths appear at each time-step. This could not be done by simply picking regions below sea-level, because this approach would, for example, cause 350 flow to "disappear" into the below-sea-level Caribou Basin of Lake Superior and/or into regions around the ice margin that were isostatically depressed below sea-level. Therefore, I developed an algorithm to define as "ocean" only those regions that were below sea-level and contiguous with the global ocean. First, I thresholded the topography-plus-ice flow-routing map, Z r , at 0, to define a binary raster map. Then, I converted this into a vector map. Using the topological analysis built into 355 GRASS GIS, I looked only for areas that touched the map edge, and chose these to be ocean. The remaining area was defined as continent and used as the flow-routing surface.

Flow routing and accumulation
Our flow-routing calculations produce maps of (1) drainage direction, (2) meteoric water discharge, (3) meltwater discharge, and (4) total water discharge. These maps are spatially-distributed, pro-360 viding 30-arcsecond-resolution time-series of the time-averaged surface water discharge that is produced by each ice-model-solid-Earth-model combination.
The flow-routing algorithm I used, r.watershed (Metz et al., 2011), uses a least-cost path search algorithm that sets it apart from other drainage analyses in that it is capable of routing water through depressions without explicitly flooding them (see, e.g., Schwanghart and Scherler, 2014 This dilemma over the Great Basin highlights a key question in how much to adjust the model by hand to fit data, versus how much to just allow the model to run. It would be possible to modify the flow-routing inputs to declare the Great Basin to be closed throughout part or all of the modeled geologic past. The approach taken here is to leave the flow-routing grid unmodified in order to 385 maintain distance between data and model (see the final paragraph in Section 2.1.3).

Distributed melt delivery to the coast
These methods also allow us to produce grids of meltwater discharge to the coast that are fullydistributed in space and time. These have been converted to GCM-input-appropriate spatial resolutions and used in preliminary work to connect ice melt and global ocean circulation and climate 390 change (Ivanović et al., 2014). This is an improvement over "hosing" experiments (e.g., Condron and Winsor, 2012), which have thus far been central to any model-based understanding of meltwater impacts on climate change. In "hosing" experiments, meltwater additions are prescribed to a patch of ocean. The fully-distributed GCM meltwater inputs create a more realistic ice-age ocean that is better-conditioned to respond to changes in meltwater inputs (that again, are distributed), 395 and can be used to connect models of ice-sheet history with their likely impacts on global climate.
Code to produce these GCM inputs is also provided in the Github repository (https://github.com/ umn-earth-surface/ice-age-rivers).   Lehner et al. (2013). b. Modern river mouths alongside the river mouth regions used to define a flow path as part of a drainage basin (Section 2.5). Names are also given in European non-English languages where they are common, alongside a non-exhaustive set of common names in local languages. Note that the mouth of the Mississippi River lies entirely outside its river mouth region. This is because the topographically-routed flow (Metz et al., 2011) follows the course of the Atchafalaya River; the present course of the Mississippi runs counter to the regional topographic gradient because it is held in its old course by the Old River Control

Drainage basins and river discharge
Structure (e.g., Harmar and Clifford, 2007 First, I create a map of regions where I expect the mouth(s) of each river system to be. Each "river 405 mouth region" (Fig. 4b) is defined by hand to accommodate the changing locations of river mouths with time as a function of sea level, ice cover, and GIA ( Figure 5. This is especially important for drainage systems such as Hudson Strait and the Saint Lawrence River, which previously played host to ice streams and now are fully or partially inundated with seawater. I then create a set of points at which each of the "major" rivers (≥1000 m 3 s −1 reconstructed 410 discharge) reaches the coast. Wherever one of these river mouths is located inside a named region ( Fig. 4), I build a drainage basin. Each drainage basin is labeled with the name of its respective river mouth region (e.g., "Mississippi" or "Columbia"). River mouth regions that contain more than one river mouth produce a named drainage basin that is an amalgamation of the drainage basins of all ≥1000 m 3 s −1 rivers in the basin; this allows us, for example, to combine the Mississippi,  (1) and (2) across each basin. Each sum, respectively, provides ice-sheet meltwater discharge (Q i ), meteoric (seasonal precipitation minus evapotranspiration) discharge (Q m ), and total river discharge (Q = Q m + Q i ). Many of these calculations provide discharges at the modern time-step that are close to the modern pre-human-impact mean river discharge 425 (Table 1, which also includes sediment discharges for reference) For some, however, the result differs somewhat due either to (1) the water balance maps not perfectly representing water delivery to the coast and/or (2) Figure 5. Sea level curves extracted from GIA models at the modern locations for the river mouths that are important to this study (Fig. 4) Kammerer (1990) a This is the sum of river inputs to Hudson Bay from pre-1975 data b Modeled following Syvitski et al. (2003) c Not including modern diversions into or out of the Great Lakes. These include the Chicago Diversion (out of the Great Lakes basin), and the Ogoki and Long Lake Diversions (into the Great Lakes basin). The drainage basin figures show the modern Great Lakes basin, but as these diversions are small compared to other sources of possible error, I have not returned the mapped drainage basin to its pre-diversion state. where there are chronological controls on ice-stream reorganization, the changes in drainage tend to produce modest and localized changes in catchment area and/or outlet location (e.g., Stokes et al., 2009;Ross et al., 2009;Ó Cofaigh et al., 2010). Better chronologies often come from direct dating of ice-marginal positions (e.g., Licciardi et al., 1999;Dyke, 2004;Gowan, 2013), fluvial deposits (e.g., Bretz, 1969;Atwater, 1984;Ridge, 1997;Benito and O'Connor, 2003;Knox, 2007;Rittenour 455 et al., 2007), and lacustrine deposits (e.g., Antevs, 1922;Kehew et al., 1994;Rayburn et al., 2005Rayburn et al., , 2007Rayburn et al., , 2011Richard et al., 2005;Breckenridge, 2007;Breckenridge et al., 2012) and these provide the fourth layer to define drainage chronologies. The fifth layer of data comes from the offshore stratigraphic record (e.g., Andrews and Tedesco, 1992;Andrews et al., 1999;Flower et al., 2004;Carlson et al., 2007;Williams et al., 2012;Maccali et al., 2013), which can be compiled into a 460 paleohydrograph that may correlate with drainage basin area (Wickert et al., 2013), but does not independently indicate where the past drainage basin margins lie.
While very detailed and well-dated geological constraints permit a quantitative approach to reconstruct past discharge, many are much more basic, simply showing when periods of high or low discharge may have existed due to changes in drainage basin area (e.g., Ridge, 1997), isotopic ratios 465 (e.g., Andrews and Dunhill, 2004), and/or proglacial lake outlets (e.g., Teller and Leverington, 2004, and later revisions). I therefore first took the simplest approach to these records, defining them as simply "high discharge" (light and dark gray on Figs. 10-15) or "low discharge" Figs. (10-15). Dark gray indicates times for which there is geological evidence for a distinct period high flow, and light gray indicates the time over which data-based inferences indicate that enhanced flow would have 470 occurred and/or when discharge would have been higher but not as high as during the time indicated by the dark gray regions. For example, Mississippi River discharge increased when Lake Agassiz flooded into it (dark gray), but was overall higher than present from the LGM until ∼12.9 ka due to ice-sheet input and a larger drainage basin area (light gray).

475
The result of the computational work is two major products: synthetic drainage basin and discharge histories. The former are summarized in four sets of maps of significant points in time since the Last Glacial Maximum (Figs. 6-9). The latter are a set of reconstructed deglacial water discharges by drainage basin (Figs. 10-15). These provide a spatially-distributed and quantified view of how rivers in North America have broadly changed over the past 20,000 years.

480
The counterpart to these computed histories are those that are driven by more direct interpretations of the data. I write, "counterpart" as opposed to "ground-truth" because these are also to some extent interpretations, based on limited data. Their key utility is that they offer a much simpler path from the raw data to the interpreted drainage basin extents and periods of high discharge, and one that can be more closely tied to geological fact. As such, disagreements between these data-derived basins and 485 those derived from models, especially where the data-driven drainage basins are tightly-constrained, most likely indicate where the models should be improved.

Comparison with geologic evidence for varying drainage basin areas and discharge
Each of the river systems in this study experienced high flow and low flow at different times. The 490 high flow could be due to meltwater inputs or drainage area expansion, and often these are linked: ice-sheet advance into mid-continent drainage basins vastly increased their area (Figs. 6-8). Figure 10 includes a brief description of each of these times, and I will expand on these here with references. The data-drive drainage histories (Section 3) through North America are updated from Wickert (2014) and are broadly consistent with the work of Licciardi et al. (1999), with radiocarbon 495 ages calibrated using IntCal13 (Reimer et al., 2013).

Drainage histories by river
The Mississippi, Hudson, and Susquehanna Rivers traded drainage inputs from the Great Lakes and their respective former ice lobes, resulting in a pattern of symmetric increases and decreases in discharge along the southern Laurentide margin. At the Last Glacial Maximum, the drainage area 500 of the Great Lakes was split between the Mississippi and Susquehanna Rivers (Ridge et al., 1991;Ridge, 1997;Licciardi et al., 1999), but by 19.9 ka, ice retreat associated with the Erie Interstade    cooling that corresponds with meltwater rerouting away from the Mississippi (Williams et al., 2012;Wickert et al., 2013) and enhanced meltwater output to the Mackenzie (Murton et al., 2010) and Saint Lawrence (Carlson et al., 2007) rivers. Model predictions here vary widely, and properly matching this continental-scale drainage rerouting away from the Mississippi basin is a feature that all successful ice models should have. The routing of flow from Lake Agassiz, at the center of the continent, was not towards the Mississippi at this time (Wickert et al., 2013), but is otherwise unclear based on current data (e.g., Carlson and Clark, 2012;Lowell et al., 2013;Breckenridge, 2015 Figure 11. Computed discharge based on the ANU ice-sheet-solid-Earth model pair (Lambeck et al., 2002), compared to the geologic record of enhanced past discharge (gray shading: dark, specific; light, inferred and/or less discharge). Notes on the geologic record can be found in Fig. 10

Rio Grande
Total discharge Meltwater discharge Meteoric water (P-ET) discharge Uncorrected modern modeled discharge

Age [ka]
Discharge [m 3 s −1 ] Figure 13. Computed discharge based on the G12/VM2 coupled ice-sheet-solid-Earth model (Gregoire et al., 2012;Peltier, 2004), compared to the geologic record of enhanced past discharge (gray shading: dark, direct record; light, inferred and/or less discharge). Notes on the geologic record can be found in Fig. 10.
Meltwater routing down the Mississippi lasted until ∼12.9 ka (Williams et al., 2012;Wickert et al., 2013), when its last major ice-sheet contributor, the Lake Agassiz drainage basin, was rerouted either 520 east to the newly-ice-free Saint Lawrence (Broecker et al., 1989;Carlson et al., 2007;Carlson and Clark, 2012) or to the Mackenzie River (Murton et al., 2010;Breckenridge, 2015). Regardless of the routing of Lake Agassiz meltwater, a significant amount of meltwater was routed to the Gulf of Saint Lawrence starting at ∼13.2 ka (Carlson et al., 2007;Rayburn et al., 2011). The Saint Lawrence continued to receive significant water inputs via lakes Agassiz and Ojibway until sometime between 525 ∼8.6 ka, the potential subglacial flood from the combined Lake Agassiz-Ojibway into Hudson Bay (Breckenridge et al., 2012), and 8.2 ka, when the ice saddle in Hudson Bay collapsed, causing the 8.2 ka event flood (Barber et al., 1999;Gregoire et al., 2012;Breckenridge et al., 2012;Stroup et al., 2013). While this was the last meltwater input to the Saint Lawrence River, meltwater from the Québec-Labrador dome flowed down the Manicouagan River and into the Gulf of Saint Lawrence 530 until ∼7.8 ka (Occhietti et al., 2004).
The record of pre-Agassiz-Ojibway flood water discharge from Hudson Bay and Hudson Strait, which was an ice-covered marine-terminating margin, comes primarily from records of ice-rafted debris in marine sediment cores. The largest of the post-LGM ice-rafted debris events was Heinrich Event 1, 18.2-17.7 ka (Rashid et al., 2003). Heinrich Event 0, corresponding to the Younger Dryas, 535 was 12.8-12.3 ka (Clark et al., 2001). Later ice-rafting events are related to ice advances in Ungava Bay; these comprise the 11.4-10.9 ka Gold Cove advance (Kaufman and Miller, 1993) and the 10.0-9.4 ka Noble Inlet advance (Metz et al., 2008). Following the ∼8.6 ka possible subglacial drainage of Lake Agassiz-Ojibway, identified by a low-water period in the varve record (Breckenridge et al., 2012), and the 8.2 ka catastrophic drainage (Barber et al., 1999), Hudson Bay was fed by a much 540 lower meltwater input from the ablating remnant Laurentide ice-caps (cf. Dyke et al., 2003;Dyke, 2004), and shrank as isostatic rebound at a time-averaged rate of 1.3 cm yr −1 (Lavoie et al., 2012) caused its shoreline to retreat.
Additional Lake McConnell outburst floods occurred at ∼13.3 (possible but uncertain) and at ∼13.1 ka. At 13.1 ka, ice retreat rerouted Glacial Lake Peace, which formerly drained into the Missouri (and hence, the Mississippi) river, towards the Mackenzie, significantly increasing its drainage area.

550
Between 13.0 and 11.5 ka, large amounts of meltwater continued to be routed down the Mackenzie River. This date range includes the likely timing of the incision of the Mackenzie River ramparts (Mackay and Mathews, 1973;Lemmen et al., 1994;Wickert, 2014), deposition of flood-transported boulders on the Mackenzie River delta (Murton et al., 2010), and a major ∼11.6 ka negative δ 18 O excursion in the Beaufort Sea (Andrews and Dunhill, 2004). After 11.5 ka, the Laurentide Ice Sheet 555 continued to retreat, feeding meltwater into the Mackenzie River system that included a 10.8-10.7 ka flood from Lake Agassiz during its Emerson phase (Fisher and Lowell, 2012). The retreating Laurentide Ice Sheet left the Mackenzie basin by ∼9.3 ka (Lemmen et al., 1994).
The Colorado River and the Rio Grande, while influenced by mountain glacier retreat, experienced different discharges largely due to the wetter climate in their drainage basins during the Last Glacial Maximum and Younger Dryas. These cool, wet periods are marked by high lake levels (e.g., Currey et al., 1990;Matsubara and Howard, 2009;McGee et al., 2012), mountain glacier advances prior to 570 the onset of the Bølling-Allerød (e.g., Owen et al., 2003;Guido et al., 2007;Orme, 2008), increased landsliding in the Rio Grande basin between ∼21.2 and ∼14.5 ka (Reneau and Dethier, 1996), and speleothem growth (Polyak et al., 2004). This pattern of climate change closely mirrors that seen in Greenland (Benson et al., 1997;Svensson et al., 2008;Asmerom et al., 2010;Wagner et al., 2010).

Comparison of data and models for river discharge 575
In order to determine how well, broadly, these different models were able to reproduce discharge, I subdivided the computed discharges into events during known high-flow periods (light and dark gray regions on Figs.10-15 and those from the low-flow remainder of the record. I first compared the mean discharges during known high-and low-flow periods (Fig. 16A). For the five ice-sheet models 31 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-8, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 February 2016 c Author(s) 2016. CC-BY 3.0 License. and nine river systems, only once did the known low-flow time produce a higher computed discharge 580 (the Hudson River for the ANU model). The ice-physics-based G12 performed the best in generating discharges during high-flow periods that were much higher than those during low-flow periods.
I then generated histograms for the high-flow and low-flow segments of the model results, and performed Kolmogorov-Smirnov tests between the pairs of distributions for each ice model and river system to test how different they are from one another (Fig. 16B). Here, the new GIA-based ICE-6G 585 performed the best and G12 performed the worst. While visual inspection of Figure 13 shows that G12 performs quite well in representing general patterns of deglacial water discharge, including the pulses of meltwater output from Hudson Strait associated that would produce the diagnostic iceberg pulses of Heinrich events (e.g., Heinrich, 1988;Andrews and MacLean, 2003;Hemming, 2004), its short (but intense) pulse down the Mackenzie River and lack of early meltwater down the Hud-590 son River reduced the average, and the averaging across the longer meltwater-production period for Hudson Strait reduced the goodness-of-fit. G12 is also affected by a too-early switch of meltwater outflow away from the Mississippi River because its climate forcing and dynamical ice response created a Québec-Labrador dome that retreats too soon.
Across all model runs, Hudson Bay was the consistent site of the worst data-model comparisons.

595
It is possible that this is due to data-collection issues (e.g., few organics) and data sparsity (few settlements and fewer GPS monuments) in this formerly fully ice-covered region, leading to a lack of good regional constraints on ice-sheet thickness near the center of the Laurentide Ice Sheet. Or, on the other hand, this could be because the models do not accurately represent the dynamics of the marine-terminating Hudson Strait ice catchment (e.g., MacAyeal, 1993).

Comparison of data and models for drainage basin extents
All of the simulated ice sheets follow the broad strokes of the expected pattern of drainage. This is as would be expected: any effort to put a realistically-shaped ice sheet onto North America will cause more flow to be routed towards the Mississippi River and less to the Saint Lawrence River.
This section covers examples of where the ice-sheet models did not match the data-driven drainage 605 basin extents, and concludes with suggestions for how to improve such fits.
At the Last Glacial Maximum, ICE-5G and G12 have the most severe misfits with data. ICE-5G sends a large fraction of Mississippi-River-bound drainage to the Hudson River, an artifact fixed in ICE-6G with its more realistic ice-dome structure. G12 at the Last Glacial Maximum has Mackenzie River drainage routed towards the Yukon River, though time-steps shortly before and after the 610 LGM have the proper drainage routing. This underscores the importance of drainage routing as a discriminating factor in ice-sheet-model quality, as it is able to detect thresholds in ice-sheet thickness that are required to maintain a realistic ice geometry. ANU and G12 produce the most realistic flow partitioning between the Susquehanna and Hudson Rivers.

32
Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-8, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 February 2016 c Author(s) 2016. CC-BY 3.0 License. Figure 16. The computed discharge time-series is segmented into values corresponding with times of known high and low discharge. a. Normalized difference between known high and low discharge times. The greatest mean normalized difference is seen in G12, the ice-physics based model of Gregoire et al. (2012), and the least in the GIA-based model of Lambeck et al. (2002). b. The Kolmogorov-Smirnov test measures the distance between one-dimensional distributions; here I use it to compute the a metric of difference between the computed discharge histograms for each river system during the generalized data-derived times of high and low discharge (both light and dark gray in Figs. 10-15). A high distance value indicates a strong separation between the range of discharges during geologically-constrained high-and low-flow periods, and hence a good match between the data and model. Here, the new GIA-based ICE-6G model performs the best, with a strong match to all but the Hudson Bay record, and the ice-physics-based G12 performs the worst, due in part to the too-early demise of its Québec-Labrador dome.

33
Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-8, 2016 Manuscript under review for journal Earth Surf. Dynam. During the early part of the Younger Dryas, the Mackenzie and Saint Lawrence Rivers both experienced periods of high discharge (de Vernal et al., 1996;Carlson et al., 2007;Murton et al., 2010). It is unclear whether their basins tapped into the center of the continent: debate continues as to whether 620 Lake Agassiz flowed towards the Mackenzie River in the northwest (e.g., Murton et al., 2010;Breckenridge, 2015), towards the Saint Lawrence River in the east (e.g., Carlson and Clark, 2012), or, as one group has proposed, became a closed basin at this time due to regional drying . Currently, the most favorable hypothesis seems to be a northwesterly route, which is consistent with the newly-remapped strandlines (Breckenridge, 2015). The drying hypothesis seems the 625 least likely, as the model simulate it  creates a closed basin when at least two parameters seem unreasonable. These are (1) evaporative loss that is generally >1.5 times today's observations, in spite of cooler temperatures and a longer lake-ice season, and (2) a temperature lapse rate that is 7.3 • C km −1 , outside of the error bars of the 5.1±1.2 • C km −1 observed on modern ice sheets and ice caps Anderson et al. (2014), and resulting in a smaller melt-producing area of 630 the Laurentide Ice Sheet. Regardless of where Lake Agassiz flowed, or whether it became a closed basin (the flow-routing presented here does not allow closed basins to exist: see Section 2.3.2), isotopic data from the Gulf of Mexico show unequivocally that meltwater discharge to the Mississippi ended ∼12.9 ka (Wickert et al., 2013). This means that the ICE-5G and ANU reconstructions miss the timing of the most major drainage rearrangement in the North American deglaciation. The G12 635 reconstruction has such intense isostatic adjustment when coupled with the VM2 solid Earth model that there is marine inundation of the entire set of Great Lakes Basins, and the Upper Missouri River in G12 and ICE-6G flows towards the north. In G12, this is due in part to the too-early retreat of the southeastern Laurentide margin. Overall, ICE-3G and ICE-6G perform the best, with the latter sending Lake Agassiz meltwater to the Mackenzie River and Beaufort Sea and water from the Laurentide 640 Great Lakes to the Gulf of Saint Lawrence, a scenario that is consistent with current geological evidence (e.g., de Vernal et al., 1996;Richard et al., 2005;Carlson et al., 2007;Murton et al., 2010;Rayburn et al., 2011;Breckenridge, 2015).
The early Holocene in northern North America was characterized by retreat of the Laurentide Ice Sheet that led to expansion of proglacial lakes, setting the stage for the short-lived and rapid 8.2 645 ka cooling event (Breckenridge et al., 2012). G12 produced a saddle collapse that led to this flood (Gregoire et al., 2012), though with a somewhat different spatial pattern than what has been mapped (Dyke et al., 2003;Dyke, 2004). None of the models perform particularly well here. ICE-3G has altogether too much ice retreat. ICE-5G, ICE-6G, and ANU have a drainage divide that splits the large east-west Saint Lawrence drainage basin, which is constrained by the dual facts that Lake Lawrence (Teller and Leverington, 2004;Breckenridge et al., 2012;Stroup et al., 2013). G12 has such a continuous drainage basin, but it flows out a gap along the southwest corner of the retreating ice-sheet, and from there, flows to the Arctic instead of to the Saint Lawrence. It is notable that ICE-5G, ICE-6G, and G12 have a tiny Mackenzie River basin due to glacial isostatic adjustment, 655 and evaluating the plausibility of this could be a goal for future field work.
5.2 Future directions: coupling with hydrologic, climate, isotopic, archaeological, and sediment transport models and observations The work presented here is a necessary starting point for future studies on the interactions between deglacial surface-water hydrology and other components of the integrated Earth system. Combining 660 these meltwater routing patterns with proglacial lakes and changing land cover will allow equations designed to predict sediment yield from large catchments to be employed for a continentalscale paleo-sediment-discharge reconstruction (following Overeem et al., 2005;Kettner and Syvitski, 2008;Pelletier, 2012;Cohen et al., 2013), with modern control provided by present-day sediment yields (Table 1). These sediment discharges can be compared with records of deposition (e.g., 665 Andrews and Dunhill, 2004;Breckenridge, 2007;Williams et al., 2010) and records of geomorphic change (e.g., Reusser et al., 2006;Knox, 2007;Anderson, 2015). The need to properly compute past lake and land cover motivates continued work with climate-and water-balance models (e.g., Collins et al., 2006;Matsubara and Howard, 2009;Liu et al., 2009;He, 2011;Blois et al., 2013).
These, together with paleogeographic reconstructions such as those presented here, can be used to 670 reconstruct areas of archaeological interest, either as changes in shoreline positions and topography  or as a wholesale landscape reconstruction that incorporates site-potential modeling (Monteleone et al., 2013;Monteleone, 2013;Dixon and Monteleone, 2014). On a global scale, several current flow-routing algorithms could be made global for better integration with ice-sheet, climate, and GIA models (Metz et al., 2011;Qin and Zhan, 2012;Braun and Bezada, 2013;Huang 675 and Lee, 2013; Schwanghart and Scherler, 2014), including the possibility to include high-resolution flow routing as part of a transient coupled model run instead of an a posteriori analysis, as is presented here. Finally, high-resolution drainage routing schemes can connect models of past climate, ice sheets, and drainage routing, to oxygen isotopes in sediment cores (Wickert et al., 2013), and the increasingly-complete collection of such records from the North American continental margin 680 (Andrews et al., 1994;de Vernal et al., 1996;Brown, 2011;Williams et al., 2012;Gibb et al., 2014;Taylor et al., 2014) is opening new possibilities in isotopic studies of whole-ice-sheet mass balance.

Conclusions
The goal of this work is to move towards a self-consistent paleogeographic framework within which models and geologic records may be quantitatively compared to build new insights into past glacial 685 35 Earth Surf. Dynam. Discuss., doi:10.5194/esurf-2016-8, 2016 Manuscript under review for journal Earth Surf. Dynam. Published: 17 February 2016 c Author(s) 2016. CC-BY 3.0 License. systems and climate-ice-sheet interactions. Most ice-sheet reconstructions are built to fit glacial isostatic adjustment measurements and moraine positions, and/or by using ice-physics-based models.
These are not calibrated to drainage basins and meltwater pathways, even though these leave behind geomorphic and stratigraphic markers that can be used to discriminate between hypothesized icethickness distributions and deglaciation histories. The dynamic deglacial drainage basins of North 690 America are significant in their own right, and the new paleohydrographs presented here can help us to better understand the ice-age legacy imprinted on the major rivers of North America. When integrated more broadly into studies of North American deglaciation, these time-varying deglacial drainage basins and paleohydrographs are poised to improve reconstructions of past climate and ice-sheet geometry.