On a map of Maine there is a big indentation midway along the coast – Penobscot Bay – the sink for the 175 km-long Penobscot River (Fig. 1). The Penobscot River heads in a series of glacially scoured lakes on an intermontane plateau, ~350 m asl (above present sea level) in Maine’s western mountains. Thence it flows through some low foothills and, near the town of Mattawaumkeag, into the Penobscot Lowland. The Penobscot River presently drains somewhat over 20 000 km² with a mean annual discharge of 400 m³/s (United States Geological Survey 2017). It is one of three large rivers reaching the Maine coast. The river’s post-glacial history in the Lowland is described in a companion paper (Hooke et al. 2016).

The last ice advance across Maine culminated ~22–20 ka (Lambeck et al. 2014). Ice then extended across the Gulf of Maine to the edge of the continental shelf. By ~16 ka, the ice margin had retreated to roughly the present coastline and, although substantial rebound had occurred, shoreline features now ~70 m asl near the present coast show that the land was still significantly depressed. Thus, with further ice retreat the sea advanced inland, inundating the Lowland, and eventually forming estuaries extending some distance into Maine’s western mountains (see marine limit in Fig. 1).

The East Branch of the Penobscot River flows in a valley which, during late glacial times, was one such estuary. The northern ~30 km of the ~150 km-long Katahdin esker (Fig. 1) runs along the western side of this valley and is rarely more than a couple of hundred metres from the river. This part of the esker was first described, albeit briefly, by Stone (1899, p. 105), who also noted “a few” expansions of the esker into “elongate plains, not very broad, but rather flat on top, sometimes enclosing kettle holes.” We have identified ten such expansions which we interpret to be deltas (Fig. 2a–d), an interpretation previously advanced by Clinch and Weddle (1989). These deltas were formed where the retreating ice paused and the subglacial esker stream delivered enough sand and gravel to the margin to form a deposit reaching sea level.

Our objectives in this paper are to describe the esker and the associated deltas in greater detail, to elucidate the rebirth of the Penobscot River as the ice retreated, and to consider the fluvial adjustments once the land emerged above the sea.

Our identifications of map units are based on topography and sedimentology. Vertical exposures being rare, the sedimentologvbnvbnvbny was commonly determined by digging shallow holes in the surface or by sampling with an 18 mm soil probe thrust to a depth of 0.1 to 0.9 m, the depth depending on the caliber and stiffness of the sediment. Sand sizes were evaluated by comparison with images on an AMSTRAT transparent sediment-size card viewed with a 10× hand lens. Gravel sizes were estimated by eye.

The primary topographic landforms are gravel ridges (the eskers) commonly exceeding 10 m in height, the deltas, and relatively planar surfaces (other than delta plains) identified either as parts of the present floodplain or as river terraces. Elevations of deltas and some terraces are central to our analysis. These were estimated from contours on 7½ minute topographic maps (see Fig. 1 for quadrangles used). Contour intervals were 3.06 m or, in the case of two northerly deltas, 6.1 m. Delta surfaces had at least 3 m of relief, and the forest cover is conducive to cartographic errors. Thus, ± 3 m is perhaps a minimum estimate of the uncertainty in these elevations.

The elevations of many terraces were also expressed in terms of their height above the level of the mean annual flood (approximately the modern floodplain). The gage height of this flood, which has a recurrence interval of ~2.33 years (Pizzuto 1986), was determined from a gaging station with a 60-year record (United States Geological Survey 2002) at Grindstone (Fig. 1; see Leopold et al. 1964, p. 63, for technique). The height of a terrace surface relative to this gage height was estimated by measuring its height above the river with a hand level and obtaining the river level on that day from the gaging station record.

In one case, heights of points on three terraces were obtained by surveying from a United States Geological Survey recoverable mark, in this case a bridge deck, with a Nikon AE-5 automatic level (Fig. 2a, AE-5 transect). The elevations of recoverable marks are given to the nearest 0.3 m, and in this case, as the bridge had been washed out, there is a further uncertainty of ~0.4 m in our estimate of where the bridge deck intersected the surface of the logging road. This results in an uncertainty of ± ~0.5 m in these elevations.

Between Medway (Fig. 1) and the coast, and extending from New Hampshire on the southwest to New Brunswick on the northeast, there are many ice-marginal deltas similar to those we mapped. Thompson et al. (1989) measured elevations of topset/foreset contacts (or in 4 cases, meltwater channels) on 69 of these. They then fit a third-degree trend surface through the elevations, obtaining a close fit to their measurements (r² = 0.98). This trend surface is an estimate of the local marine limit as a function of distance from the coast and time. While the deltas were being built ice was retreating, the crust was rebounding, and global mean sea level was rising, the latter two at rates of a metre or two per hundred years. Thus, the trend surface does not represent sea level at any specific time in the past. In order to extrapolate this trend surface north of Medway, we fit a parabola to it along a generally south-north line roughly coincident with the Penobscot River.

Samples of sand for optically stimulated luminescence (OSL) dating were collected from the sides of small gravel pits in two terraces along the AE-5 transect (Fig. 2a). We first dug back several decimetres to expose undisturbed stratigraphy and then hammered 2” × 7” PVC tubes into the resulting face. To avoid recent soil disturbance, samples were collected from depths >1.5 m. After removing a tube from the face, excess sand was trimmed leaving a sand surface flush with the tube ends to minimize disturbance during shipping. The ends were then covered with multiple layers of duct tape. The samples were processed at the Luminescence Geochronology Laboratory at the University of Nebraska-Lincoln. Equivalent dose (D_e) values were obtained using the single aliquot regenerative (SAR) method (Murray and Wintle 2000) on 90–150 μm quartz grains. OSL measurements were conducted on a Risø model DA-20 luminescence reader, and a preheat plateau test (Wintle and Murray 2006) was used to determine the optimum preheat temperature. Aliquots were rejected if their recycling ratios were >10% or if they had measurable signals during stimulation with IR diodes. A minimum of 30 accepted aliquots were generated for each sample, and final age estimates were calculated using the Central Age Model of Galbraith et al. (1999). Environmental dose rate estimates were based on elemental concentrations of bulk sediment taken from within ~30 cm of the OSL sample. These bulk sediment samples were analyzed for concentrations of K, U, and Th using a high-resolution gamma spectrometer. The cosmogenic component of the dose rate values was calculated using equations from Prescott and Hutton (1994). The long-term moisture content was estimated to be 8%. The final dose rates were calculated following equations froM Aitken (1998). Age estimates, dose rate values, and equivalent dose values are shown in Table 1.

We constructed a time-distance plot of ice sheet retreat (Fig. 3) along a transect running from the coast, through Ellsworth, and ~10 km west of Lower Togue Pond (Fig. 1). The transect is roughly normal to contours of deglaciation age drawn by Borns et al. (2004). All but two of the dates shown in Figure 3 are ¹⁴C dates reported by Borns et al. (2004) and calibrated using CAILB 7.0 (Stuiver et al. 2014). The other two are from this study. Most of the dates are from sites lying within ~50 km of the transect and projected onto the line of transect along lines parallel to the contours of deglaciation age. Much of the scatter is likely due to errors in interpreting the context of the dates and in projecting the sites onto the transect. Retreat rates are based on the eight points shown with filled symbols and listed in Table 2.

Several of the dates were on marine organisms and required a marine reservoir correction. Reservoir corrections vary both spatially and temporally. Along the Atlantic coast from New Jersey to Newfoundland, McNeely et al. (2006) report values ranging from 430 to 730 years. By comparing dates on shells and on associated logs in the Presumpscot Formation near Portland, Maine, Thompson et al. (2011) and found a reservoir age of 1000 years. Following Borns et al. (2004) and consistent with a rough estimate of the median value obtained by McNeely et al. (2006) for the Atlantic coast between New Jersey and Newfoundland, we used 600 years.

Hooke et al. (2016) suggested that two prominent moraines near the coast, the Pineo Ridge (Fig. 1) and nearby Pond Ridge moraines, were likely formed during minor cold events at 15.6 and 16.1 ka, respectively. These events were recorded as oxygen isotope minima in the NGRIP (North Greenland Ice Core Project) ice core (Fig. 4) (GICC05 chronology from Andersen et al. 2006). Such cold events would likely have resulted in pauses in ice retreat, or possibly slight readvances. We thus infer that they were responsible for these moraines.

The dates of the cold events are generally consistent with available¹⁴C dates constraining the ages of these moraines (Table 2). At Dennison Point, 3 km outside the Pond Ridge moraine, a Macoma calcarea shell closely associated with ice-proximal grounding-line deposits was dated to 16.1 ± 0.1 ka, the same as our estimate of the age of the Pond Ridge moraine. On the distal side of the Pond Ridge moraine at Turner Brook a Nucula sp. shell from glaciomarine sediments interbedded with coarse ice-proximal units was dated to 15.9 ± 0.1 ka. Finally, on the proximal side of the Pond Ridge moraine at Sprague Neck a Portlandia arctica shell in glaciomarine facies was dated to 15.2 ± 0.1 ka, yielding a minimum age for the moraine. Similarly, our age for the Pineo Ridge moraine (15.6 ka) is consistent with basal¹⁴C dates from Marks and Lily lakes, a few kilometres inside the moraine, that yield minimum ages for the moraine of 15.1 ± 0.1 and 15.2 ± 0.1 ka, respectively.

Further north, near Lincoln (Fig. 1) a Nucula sp. shell from near the base of some glaciomarine mud yielded an age of 15.1 ka (T2 R8 in Fig. 3 and Table 2) which Dorion et al. (2001, p. 235) interpreted to be a “close minimum” age for deglaciation. This suggests a retreat from the Pineo Ridge moraine to this point, a distance of ~100 ± 10 km, in ~500 ± 140 years, or 200 ± 60 m/yr. This is somewhat faster than, but not inconsistent with a mean rate (~150 m/yr) deduced by Eusden (2014) based on a study of the spacing of a plethora of likely roughly annual minor moraines in the southern Penobscot Lowland.

Still further north, Davis and Davis (1980) obtained a date of 13.5 ka on terrestrial materials 0.5 m above the base of a core from Lower Togue Pond (Figs. 1, 3). Extrapolating to the base of the core, assuming a uniform sedimentation rate over the 13 m length of the core, suggests deglaciation at 14.1 ka. Two other nearly identical dates (Anderson et al. 1986) on plant macrofossils from 3.3 and 4.3 m above the base of a core from Upper South Branch Pond (Figs. 1, 3) yield a mean age of 12.9 ± 0.2 ka and a deglaciation age, based on crude curve fitting to scattered dates, of 13.9 ka. Thus, the Seboeis delta was likely formed at ~14 ka. The retreat rate from Lincoln to the Seboeis delta, a distance of 65 km, was then ~76 m/yr (Fig. 3).

The rate of change of relative sea level at any given location on the globe is a function of three processes: (1) Glacioisostatic rebound, (2) changes in the gravitational attraction exerted by the ice sheet on the sea, and (3) ice-sheet decay. Rebound is expected to be essentially continuous with only minor fluctuations; the rate decreases through time, asymptotically approaching zero. The gravitational attraction exerted on the sea by the ice sheet changes with changes in the mass of the ice sheet and its distance from the sea (Clark 1976). Ice sheet decay is driven by climate and is irregular. Relative sea level decreases when rebound plus the gravitational effect exceed meltwater addition to the ocean. Pauses in such decreases reflect times when meltwater addition accelerates to the point where it equals the other two.

At about the time of the rebirth of the Penobscot River, sea level was rising unusually rapidly in what Fairbanks (1989) has termed meltwater Pulse 1A (MWP-1A). Lam-beck et al. (2014 fig. S4B) date this to 14.5–14.1 ka (Fig. 4) but other sources suggest that the peak may have been later—at ~13.8 ka (Gornitz 2009). About 3 k.y. later there seems to have been a second pulse—MWP-1B (Fairbanks 1989). The rate of rise during these pulses varied spatially, but was typically a few tens of millimetres per year.

If MWP-1A occurred during retreat of the ice margin up the East Branch, as our dating suggests, the consequent rise in sea level would result in delta elevations above our extrapolation of Thompson et al.’s (1989) data (Fig. 5). However, if the rebirth of the Penobscot River occurred a few centuries earlier than 14 ka, or MWP-1A was somewhat later, the pauses in the fall of relative sea level represented by the Deasey terrace, and perhaps by the various shoreline nips, could be attributable to it. Then the Owen terrace could be related to MWP-1B. As the accuracy of dating techniques improves in the future, it may become possible to more definitively associate such geomorphic features with sea level history.

The reborn Penobscot River must have extended itself rapidly, both headward and seaward. Although the Bølling warm peak had passed by 14 ka (Fig. 4) and the climate was descending into the Younger Dryas, it is safe to assume that the ice sheet had not fully adjusted to the rapid warming of the previous millennium and was still retreating. Concurrently, the river was being extended southwards as isostatic rebound and relaxation of the gravitational pull of the ice sheet on the sea generally outpaced global sea level rise from melting ice sheets, and relative sea level was falling (Kelley et al. 2013).

Further south along the East Branch, there are several terraces and discontinuous subhorizontal gravel-mantled surfaces that may be parts of terraces. Some of these are shown on earlier reconnaissance maps (Newman 1986; Newman and Genes 1986). The surfaces were not mapped in detail and are not shown on maps herein. They lie, 1 to as much as 12 m above the modern floodplain, and are underlain by gravel consisting of rounded pebbles and small cobbles. In some cases, the gravel is overlain by 1.5–2 m of fining-upward silty fine to medium sand.

The highest such surface that we found is ~4 km north of Medway on the east side of the river. Exposures in the track of a woods road suggest that it is composed of gravelly fine sand. It is ~20 m below the marine limit. Its origin, and that of two other gravel-mantled surfaces at similar heights shown on Newman’s (1986) map, are problematic. They do not appear to be extensive enough to represent a valley fill associated with the rebirth of the river, but they are well below the marine limit, and thus unlikely to be ice-marginal subaerial kames. Lower terraces, those 1 to 2 m above the present floodplain, are likely straths cut during Holocene downcutting.

The modern floodplain is inset into these lower terraces. It, too, is discontinuous and is underlain by alternating beds of pebble gravel, fine to medium sand, and poorly sorted silt and very fine sand, commonly fining upward. At a depth of ~2 m there is an abrupt change to pebble-small cobble braidplain gravel. This gravel bed appears to extend all the way to the present coast (Hooke et al. 2016). The gravel was likely eroded from the eastern slopes of the esker and the deltas in early post-glacial time when discharges were higher and while the Penobscot River was extending itself southward.

As the Laurentide ice sheet retreated across the Penobscot Lowland between 15.7 and 15.1 ka and up the East Branch between 15.1 and 14 ka, the area was isostatically depressed and inundated by the sea. Retreat across the Lowland was rapid, averaging ~200 m/a; retreat up the East Branch slowed to ~76 m/a. During pauses in the retreat up the East Branch ten deltas were formed. As the deltas occupy a 32-km stretch of the valley, they were likely deposited over a time span of ~420 years (32 km/76 m·a^-1).

At the northing of the low Seboeis delta, land became exposed in the bottom of the valley of the East Branch at ~14 ka: the Penobscot River was reborn. Between ~14 and ~12 ka, the Penobscot extended itself, both to the north as the ice retreated, and to the south as relative sea level fell, building a braidplain in the process—the 118-m braidplain. Once relative sea level had fallen ~4 m, it appears to have stabilized briefly, resulting in regrading of the 118-m surface to form the Deasey surface and possibly in erosion of the shoreline nips in several deltas.

Discharges were likely larger than at present during this time span as ice was still melting in the headwaters and the diversion of flow from the Penobscot basin into the Kennebec River at Moosehead Lake had not yet occurred. Differential rebound was incomplete so the river gradient was lower. Despite this lower gradient, the high discharges transported gravel eroded from the eastern slopes of the esker and deltas and deposited it on strath surfaces all the way to the coast.

After a further drop of ~3 m in relative sea level, another brief period of stability led to regrading that formed the Owen surface at ~112 m, leaving the Deasey surface as a terrace.

Although the periods of stable relative sea-level that resulted in the remarkably uniform heights of the deltas of the East Branch and in the Deasey and Owen terraces apparently occurred close to the time of MWP-1A and 1B, it is presently impossible to relate these periods of stability to specific events in the global sea level curve.

Initial downcutting from the Owen surface was likely driven by falling sea level, but more recently a decline in sediment load and an increase in the river gradient due to isostatic tilting have probably played a role. Lowered discharge at ~10 ka resulted in a transition to fine-grained over-bank sedimentation.