Glacial isostatic adjustment deflects the path of the ancestral Hudson River

The Laurentide Ice Sheet reached its maximum extent over the last glacial cycle at the Last Glacial Maximum (LGM, 26 ka). Evidence of ice-sheet extent during the glaciation phase preceding the LGM is obscured by the subsequent retreat of the ice sheet, and is poorly preserved in the geological record (Clark et al., 1993; Dyke et al., 2002; Kleman et al., 2010). Geochronological constraints on ice-sheet history are difficult to obtain because much of the last ice age is beyond the age limit of radiocarbon dating, and the remaining geologic observations are characterized by sparse dates, using methods such as optical stimulated luminescence (OSL) or cosmogenic nuclides (Briner et al., 2006; Dalton et al., 2016). These limitations extend to other regions that were glaciated across the last ice-age cycle, and they pose a serious challenge to studies of glacial isostatic adjustment (GIA) prior to the LGM, in addition to hindering our understanding of the dynamics of ice-sheet growth. However, it should be possible to gather insight into ice-sheet histories by examining well-dated evidence of landscape evolution in unglaciated regions proximal to continental glaciation. GIA produces geographically variable crustal uplift and subsidence rates that reach many tens of meters per thousand years, a level sufficient to control river drainage patterns (Wallinga et al., 2004; Wickert et al., 2013; Wickert 2016) and delta morphology (Whitehouse et al., 2007).

The diversion of the Hudson River (northeastern North America) during the latest phase of Marine Isotope Stage 3 (MIS 3; 60–26 ka; Siddall et al. 2008) is an example of landscape evolution that has been connected—albeit speculatively—to crustal deformation related to the growth of the Laurentide Ice Sheet and, in particular, the dynamics of its peripheral bulge (Knebel et al., 1979; Carey et al., 2005). The river was diverted eastward from its ancestral channel (dashed black line in Fig. 1), which has been dated to 40–30 ka using radiocarbon (¹⁴C) and amino acid racemization dating (Knebel et al., 1979; Sheridan et al, 2000; Carey et al., 2005) and mapped by seismic reflection surveys (Carey et al., 2005; Christensen et al., 2013; Santra et al., 2013). Though this river diversion at ca. 30 ka has been attributed to ice-age processes, this hypothesis has not been tested using process-based models of either GIA or landscape evolution.

Predicting how GIA modulates these river processes requires knowledge of the history of ice growth over North America prior to the LGM. Glacial geologists have mapped out terminal moraines tracing the maximum extent of the Laurentide Ice Sheet at the LGM: the Erie and Lake Champlain lobes entered Pennsylvania from the northwest and northeast, respectively, forming a triangular-shaped ice-free region between the two ice lobes (Braun 2004). Ice advanced to its maximum position in southern New England at 27–20 ka (Balco and Schaefer 2006); however, little evidence remains in New England of ice extent during the build-up phase (Corbett et al., 2017). Pre-LGM configurations of the Laurentide Ice Sheet have been constructed by Clark et al. (1993) on the basis of extensive stratigraphic records, and by Kleman et al. (2010) using a wide range of geomorphic evidence (eskers, ribbed moraines, ice-flow traces, and glacial striations). Nevertheless, given the uncertainties in chronology discussed above, it has not been possible to identify a spatially and temporally precise Laurentide Ice Sheet configuration during MIS 3. Recently, Dalton et al. (2016) inferred delayed glaciation of the eastern sector of the Laurentide Ice Sheet by dating non-glacial deposits in the region to 50–37 ka using OSL, uranium-thorium (U-Th), and ¹⁴C chronometers. They concluded that the Hudson Bay Lowlands was possibly ice-free from MIS 5a (100 ka) through mid-MIS 3 (ca. 40 ka).

In this study, we predict crustal deformation using GIA simulations based on a suite of possible ice-sheet histories and viscoelastic Earth structures, and use the predicted changes in topography to drive a landscape evolution model that allows river networks to evolve under spatially variable uplift. We use these results to explore the extent to which GIA may have controlled the diversion of the Hudson River.

Ice-sheet growth and decay drive a complex pattern of sea-level (or equivalently, topographic) change. In our simulations, we performed sea-level calculations based on the theory and pseudo-spectral algorithm described by Kendall et al. (2005) with a spherical harmonic truncation at degree and order 256. In our primary GIA calculations, we adopt an Earth model with a lithospheric thickness of 96 km, and upper and lower mantle viscosities of 0.5 × 10²¹ Pa·s and 1.5 × 10²² Pa·s, respectively. This model is consistent with recent GIA studies involving data from the United States Atlantic coastal plain (Potter and Lambeck 2003; Creveling et al., 2017; Pico et al., 2017). In the GSA Data Repository¹, we consider the sensitivity of our results to these Earth model parameters.

We constructed a series of ice histories spanning the last glacial cycle that are used as input to the GIA simulations. In one case, we adopt the global ice history ICE-5G (Peltier and Fairbanks, 2006) across the last glacial cycle with two important modifications. The original Version 1.2 of the ICE-5G history assumes that prior to 26 ka, the ice sheet perimeter was fixed at the maximum perimeter. Our version adopts the same global mean sea-level variation as ICE-5G (Fig. 2A; dotted black line), but we make the common assumption that the geographic distribution of ice before 26 ka is the same as that afterwards, during the last deglaciation, for a given global ice volume (Raymo et al., 2011). Furthermore, we revise the ICE-5G LGM ice distribution to spatially resolve ice lobes described in the Dyke et al. (2002) reconstruction (Briner et al., 2006; Fig. DR1A in the Data Repository). We refer to this model as “modified ICE-5G”.

Many previous inferences of global mean sea level during the last glacial phase, and particularly during MIS 3, reconstruct higher sea levels (i.e., smaller global ice volume) than the ICE-5G history (e.g., Lambeck and Chappell, 2001; Pico et al., 2016; Creveling et al., 2017). Moreover, a recent GIA analysis demonstrated that anomalously high sea-level markers dated to 50–35 ka along the U.S. mid-Atlantic coast are reconciled by adopting a significantly reduced eastern Laurentide Ice Sheet, and concluded that this sector underwent especially rapid growth from 44 to 26 ka (Pico et al., 2017). Motivated by these studies, as well as non-glacial deposits dated to mid-MIS 3 in eastern Laurentia (Dalton et al., 2016), we constructed an alternate ice history, which we denote as ICE-PC2. This ice history differs from “modified ICE-5G” by adopting a global mean sea-level history consistent with recent constraints on peak sea-level highstands across the glacial phase. In particular, these constraints suggest that global ice volumes reached above −20 m global mean sea level during both MIS 5a and 5c (ca. 100 ka and ca. 80 ka, respectively; Creveling et al., 2017) and that they tripled from 44 ka to 26 ka (Pico et al., 2016). The resulting global mean sea-level curve is shown in Figure 2A. In ICE-PC2, the eastern sector of the model Laurentide Ice Sheet is ice-free from 80 to 44 ka, and this region glaciates rapidly leading into the LGM from 44 ka to 26 ka (Fig. DR1D).

Finally, we constructed three other ice models that provide variations in the ICE-PC2 history to explore the sensitivity of river network evolution to the modeled ice history. These are denoted as ICE-PC, ICE-PC3, and ICE-PC4, and they are described in detail in the Data Repository.

To begin, we predicted the change in topography over the last glacial cycle using the “modified ICE-5G” model (black dotted line in Fig. 2A; Peltier and Fairbanks, 2006). The global mean sea-level curve associated with the modified ICE-5G history is characterized by ice volumes that reach near-LGM values by ca. 65 ka; i.e., prior to the beginning of MIS 3 (Fig. 2A). Because the solid Earth response to pre-65 ka ice-load variations has largely equilibrated by 30 ka, our modified ICE-5G–based calculation predicts low rates of crustal deformation (∼1 m/k.y.) and no pronounced regional west-to-east tilting from 32 ka to 26 ka (Fig. 2B). The observed eastward diversion of the Hudson River during this time period implies contemporaneous deformation in the region, with a decreasing gradient in crustal uplift rates from west to east, a pattern that is not predicted with the modified ICE-5G history.

We next performed GIA simulations adopting the ice history ICE-PC2, characterized by a rapid drop in global mean sea level from 44 to 26 ka (blue line in Fig. 2A). The average rate of change in topography predicted by this simulation from 32 to 26 ka is characterized by west-to-east tilting of the crust, which is particularly pronounced in the region bounding the predicted shorelines at 26 and 32 ka (dotted lines in Fig. 2C; Fig. DR2).

Next, we explored how the landscape responds to the predicted GIA-induced deformation fields. Rivers respond through increased erosion or deposition to changes in channel slope at the shoreline, and numerous studies have focused on the response of coastal landscapes to globally uniform sea-level variations (Meijer 2002; Fagherazzi et al., 2004). Here, however, we focus on the response of rivers to GIA-induced crustal deformation. To highlight these effects, our landscape evolution simulations drive river network evolution only through fluvial processes in both transport- and detachment-limited rivers, and do not include coastal processes such as delta formation.

We performed simulations of landscape evolution using a generalized initial topography in the Hudson River region at the onset of the simulation (32 ka) because topography at this time was dissimilar to that of today. In particular, at ca. 13 ka, the Hudson River drained a catastrophic flood from pro-glacial Lake Iroquois, which led to up to 40 m of incision in the present-day channel, and the deposition of debris flows in the surrounding regions of the shelf (Uchupi et al., 2001; Rayburn et al., 2005). As a consequence, maps of present-day elevation are a poor representation of the landscape at the time of the diversion of the Hudson River. Our model is instead initiated using a regional synthetic reconstruction of topography on a grid with a uniform spacing of 30 arc-seconds (∼1 km; Fig. 3A; see the Data Repository) that is characterized by a north-south gradient similar to the general trend on the Hudson shelf (1 m/km); this yields paleo-Hudson River drainage that flows from north to south, consistent with paleoflow directions at ca. 30 ka recorded in the stratigraphic record (i.e., Carey et al., 2005; Fig. 1). Our goal in establishing the initial field of topography described above is to assess the general trends in the evolution of a north-to-south–flowing river in this location when subject to deformation due to GIA.

The location of river paths on the synthetic topography field is shown in Figure 3B, and the largest rivers (measured by drainage area > 10⁸ m²) are shown in green. There are nine such rivers and ten “intermediate” (> 5 × 10⁷ m²) rivers (shown in blue in Fig. 3B). To explore whether GIA may have been responsible for the diversion of the Hudson River and whether this connection favors certain ice histories, we focus on the response of this synthetic river network to topographic changes predicted by each of the GIA simulations. Our measure of performance in this regard is the number of large and intermediate rivers diverted during the period 32 ka to 26 ka.

To quantify river network responses to GIA-induced crustal deformation, we adopt the Channel-Hillslope Integrated Landscape Development (CHILD) landscape evolution model, which includes parameterizations for both detachment-limited and transport-limited processes (Tucker et al. 2001; Tucker, 2010), both of which may influence river evolution in the study region. We drove CHILD simulations by applying time-variable uplift rates drawn from the GIA predictions from 32 to 26 ka (Fig. DR2) to the initial topography, and computed the responses in the resulting river network. The parameter values adopted for these simulations are listed in Table DR1, and we review the fluvial and hillslope processes applied in CHILD and describe additional simulations that vary parameter choices in the Data Repository. These sensitivity tests revealed that the number of modeled river diversions is relatively insensitive to the relative dominance of detachment-limited or transport-limited conditions.

Driving the landscape evolution model with the vertical displacement field produced by ice history ICE-PC2 (Fig. 2C), and adopting the erodibility parameters described in the Data Repository, predicts that all nine large rivers and six intermediate rivers are diverted eastward by 6 k.y. after the onset of the simulation at 32 ka (a total of 19 sites of eastward diversion and 2 sites of westward diversion along 17 rivers; gray circles in Fig. 4B). This trend is consistent with the eastward diversion of the Hudson River preserved in the geological record. In contrast, only six large rivers and one intermediate river are diverted eastward when the simulation is driven by deformation associated with the modified ICE-5G simulation using the same Earth model and erodibility parameters (a total of nine sites of eastward diversion and one site of westward diversion along 16 rivers; gray circles in Fig. 4A).

Analogous predictions for additional ice histories and Earth models are shown in Figures DR3–DR6. We note that the predicted west-to-east pattern of deformation does not reflect the broader geometry of the Laurentide Ice Sheet peripheral bulge, which spans the greater part of the United States, but is due to a superimposed, smaller-scale peripheral bulge formed by the Erie ice lobe (Fig. DR1A; Dyke et al., 2002). The tilting caused by the late and rapid growth of the Laurentide Ice Sheet and its associated ice lobes could have been sufficient to have affected other rivers along the U.S. east coast, including the Delaware River, a connection that warrants further study.

Previous studies have speculated that the eastward diversion of the Hudson River at ca. 30 ka was driven by changes in topography associated with ongoing GIA (e.g., Knebel et al., 1979). Our modeling of both GIA and the associated river network evolution provides strong quantitative support for this hypothesis. An eastward diversion requires that the GIA field introduce a west-to-east tilting of topography. Our simulations indicate that this tilting would have been a natural consequence of a rapid, late growth of the eastern Laurentide Ice Sheet before the LGM. These simulations suggest that the Hudson River likely was diverted during the sea-level fall as the shoreline passed through the region of maximum west-east tilting from 32 to 26 ka. Our simulations ignore shallow coastal processes, including backwater effects and autogenic avulsions; incorporating these effects may increase the likelihood of diversions.

Our results demonstrate that solid-Earth deformation over Pleistocene glacial cycles can be an important influence on landscape evolution in regions without glacial cover, particularly at the periphery of former ice sheets, where information about past river evolution can refine estimates of glaciation rates.

We appreciate insightful comments from Andrew Wickert, Torbjorn Tornqvist, and anonymous reviewers. We acknowledge support from the National Science Foundation (NSF), Graduate Research Fellowships Program (Pico), NSF grant EAR-1527351 (Mitrovica), Harvard University (Mitrovica), NSF grant EAR-1525922 (Ferrier), and American Chemical Society Petroleum Research Fund 58209-DNI8 (Ferrier).