Atmospherically produced beryllium-10 in annually laminated late-glacial sediments of the North American Varve Chronology

We attempt to synchronize the North American Varve Chronology (NAVC) with ice core and calendar year timescales by comparing records of atmospherically produced ¹⁰Be fallout in the NAVC and in ice cores. The North American Varve Chronology (NAVC) is a sequence of 5659 varves deposited in a series of proglacial lakes adjacent to the southeast margin of the retreating Laurentide Ice Sheet between approximately 18 200 and 12 500 years before present. Because properties of NAVC varves are related to climate, the NAVC is also a climate proxy record with annual resolution, and our overall goal is to place the NAVC and ice core records on the same timescale to facilitate high-resolution correlation of climate proxy variations in both. Total ¹⁰Be concentrations in NAVC sediments are within the range of those observed in other lacustrine records of ¹⁰Be fallout, but ⁹Be and ¹⁰Be concentrations considered together show that the majority of ¹⁰Be is present in glacial sediment when it enters the lake, and only a minority of total ¹⁰Be derives from atmospheric fallout at the time of sediment deposition. Because of this, an initial experiment to determine whether or not ¹⁰Be fallout variations were recorded in NAVC sediments by attempting to observe the characteristic 11-year solar cycle in short varve sections sampled at high resolution was inconclusive: short-period variations at the expected magnitude of this cycle were not distinguishable from measurement scatter. On the other hand, longer varve sequences sampled at decadal resolution display centennial-period variations in reconstructed ¹⁰Be fallout that have similar properties as coeval ¹⁰Be fallout variations recorded in ice core records. These are most prominent in glacial sections of the NAVC that were deposited in proglacial lakes and are suppressed in paraglacial sections of the NAVC that were deposited in lakes lacking direct glacial sediment input. We attribute this difference to the fact that buffering of ¹⁰Be fallout by soil adsorption can filter out short-period variations in an entirely deglaciated watershed, but such buffering cannot occur in the ablation zone of an ice sheet. This implies that proglacial lakes whose watershed is mostly glacial may effectively record ¹⁰Be fallout variations. We attempted to match centennial-period variations in reconstructed ¹⁰Be fallout flux from two segments of the NAVC with ice core fallout records. For both records, it is possible to obtain matches that result in acceptable correlation between NAVC and ice core ¹⁰Be fallout records, but the best-fitting matches for the two segments disagree, and only one of them is consistent with independent calendar year calibrations of the NAVC and therefore potentially valid. This leaves several remaining ambiguities in whether or not ¹⁰Be fallout variations can, in fact, be used for synchronizing NAVC and ice core timescales, but these could most likely be resolved by higher-resolution and replicate ¹⁰Be measurements on targeted sections of the NAVC.

1 Introduction

We describe measurements of atmospherically produced ¹⁰Be in annually laminated sediments of the North American Varve Chronology (NAVC; Ridge et al., 2012). These sediments were deposited in an initially proglacial and subsequently paraglacial lake adjacent to the southeast margin of the retreating Laurentide Ice Sheet (LIS) between 18 200 and 12 500 years BP. The NAVC records events taking place at the ice sheet margin during this time, including the position and retreat rate of the ice margin, meltwater and sediment production, and proglacial lake outburst floods (Ridge, 2004, 2012; Ridge et al., 2012). Because some of these events are related to climate (Ridge et al., 2012; Thompson et al., 2017), the NAVC is also a climate proxy record with annual resolution.

Our overall goal in investigating ¹⁰Be deposition in NAVC sediments is to synchronize this record with Greenland ice core records that provide a primary template for our understanding of Northern Hemisphere climate change during the last deglaciation. The NAVC is a floating varve chronology that is anchored to the absolute timescale by radiocarbon dating of plant macrofossils deposited in individual varves, and the uncertainty in this calibration is estimated to be no better than ±200 years (see Sect. 1.2). When changes in ice margin positions and meltwater flux recorded in the NAVC are compared with climate records in Greenland ice cores, it is evident that many events recorded by the NAVC appear to correlate with rapid temperature changes in Greenland, for example those occurring at the boundaries of Greenland stadials and interstadials GS2a, GI-1e through 1a, and GS-1 (Lowe et al., 2008). Improved precision in correlating these two records could therefore potentially illuminate important aspects of ice sheet–climate feedbacks during rapid climate changes that took place during deglaciation, for example the phasing relation between atmospheric temperature change in Greenland and the retreat rate of the southern margin of the Laurentide Ice Sheet in New England or the effect of proglacial lake drainage into the North Atlantic Ocean on regional surface climate.

In principle, it should be possible to synchronize these two climate records by comparing deposition rates of atmospherically produced ¹⁰Be. ¹⁰Be is produced in the atmosphere by cosmic-ray bombardment of target nuclei, primarily N and O, and delivered to the surface by precipitation or dry fallout (Lal and Peters, 1967; Lal, 1987). The ¹⁰Be production rate varies significantly with changes in solar activity and geomagnetic field properties, which are global in nature. Although the deposition rate is further mediated to some extent by regional atmospheric circulation, variations in the ¹⁰Be fallout rate are, in general, globally synchronous. Records of the ¹⁰Be deposition rate in ice cores and sediments have been extensively used to reconstruct past solar and geomagnetic changes (e.g., Bard et al., 1997; Muscheler et al., 2007; Steinhilber et al., 2012) as well as to synchronize different records (e.g., Adolphi et al., 2018, and references therein). Beryllium-10 deposition rates have been measured in ice cores in Greenland and elsewhere in many studies (these are reviewed and compiled in Adolphi et al., 2018) and display significant centennial-scale variability during the time period represented by the NAVC.

Our goal here is to determine whether it is possible to generate a similar record of ¹⁰Be deposition from NAVC sediments that can be correlated with the Greenland records, thereby potentially improving the absolute dating of events recorded in the NAVC and our ability to relate them to climate changes recorded in ice cores. To pursue this, we describe (i) basic observations relating to ¹⁰Be deposition and systematics in glacial and paraglacial lake sediments of the NAVC, (ii) experiments designed to determine whether or not variations in ¹⁰Be fallout due to solar variability are recorded in the NAVC, and (iii) a comparison between ice core ¹⁰Be fallout records and ¹⁰Be fallout records from portions of the NAVC.

1.1 The North American Varve Chronology

The NAVC consists of annually laminated sediments that were deposited in several lakes adjacent to the margin of the retreating Laurentide Ice Sheet (LIS) in New York and New England between approximately 18 200 and 12 500 years ago (Figs. 1, 2). These sediments occur in numerous stratigraphic sections throughout the region that were originally correlated and assembled into several floating varve chronologies by Antevs (1922, 1928). These have subsequently been consolidated, corrected in a few places, and extended into a single 5659-year sequence by Ridge et al. (2012). In this section we briefly summarize the stratigraphy and sedimentology of NAVC sediments as they are relevant to this study; a complete review can be found in Ridge et al. (2012) and Ridge (2012).

The majority of the NAVC is composed of sediments deposited in proglacial lakes that occupied the north–south-trending Connecticut River Valley from central Connecticut to northern Vermont as the margin of the LIS receded northward. Although the extent of Connecticut Valley lakes varied with ice margin position, glacioisostatic rebound, and changes in outlet position and lake level, some portion of the lake system was continuously present and in contact with the ice margin until approximately 13 700 years ago, at which point changes in the sedimentological characteristics of the varves show that the lake was no longer receiving direct glacial meltwater. The complete chronology includes ∼ 4200 glacial varves overlain by ∼ 1500 paraglacial varves; the glacial–paraglacial transition is not isochronous across the entire lake system and is also gradational over ∼ 100–200 years at some sites.

Figure 1 Map of New England and adjacent parts of the eastern United States and Canada showing the location and extent of the NAVC (Ridge et al., 2012). The NAVC consists of sediments from several proglacial lakes that existed at different times. The dark shaded areas show the maximum extent of major proglacial lakes present during deglaciation, although lake boundaries were time-transgressive throughout deglaciation and these do not represent exact lake boundaries at a specific time. Labeled triangles show the location of sections sampled for this study (the abbreviations follow Table 1), with the exception of 930PN, which is located off the map to the west in central New York. Nearly the entire time span of the NAVC is represented by sediments from Glacial Lake Hitchcock, and sequences from other nearby glacial lakes are correlated with the Lake Hitchcock sequence as shown in the stratigraphic diagram at right. Although lake sediments at any particular point span a range of NAVC years, deglaciation of the area proceeded from south to north, so, in general, younger (higher-numbered) varves are found at more northerly sites (Fig. 2).

Table 1 Characteristics and purpose of varves and sets of varves sampled for Be measurement.

^* Varves at 930PN are unusual in that they are derived from a limestone source and predominantly carbonate.

Glacial varves in the NAVC are, in general, substantially thicker than paraglacial varves and are derived predominantly from direct glacial meltwater and sediment production, although paraglacial sources including surface runoff from the recently deglaciated landscape adjacent to the lake may contribute. At any one location in the lake, the contribution of paraglacial sediment to glacial varves increases upsection as the distance to the receding ice margin increases. Glacial varves, especially in areas proximal to the glacier, have melt season (summer) layers composed of a stack of micrograded beds that often fine upwards. Summer layers grade into a non-melt season (winter) layer composed of nearly pure clay, 60 % or more of which is typically less than 1 µm in size. Glacial varve sections also occasionally have extremely thick varve couplets (over 50 cm) with single thick graded beds produced by the sudden release of water from proglacial lakes impounded in tributary valleys (see, for example, Thompson et al., 2017). Paraglacial varves are thinner, typically millimeter scale, and composed predominantly of glacial till, outwash, and other ice-marginal sediments that were deposited on the surrounding landscape and subsequently eroded and transported to the lake by surface runoff after deglaciation.

Paraglacial varves in the NAVC are considered to be varves formed with zero or negligible direct glacial input, although the transition between glacial and paraglacial varves is typically gradational. They are composed of silty and fine sandy summer layers, and may have occasional graded fine sand beds much thicker than a typical summer layer that represent unusual runoff or flood events that washed sandy sediment onto the lake floor. Winter layers in paraglacial varves have diffuse contacts with summer layers below and sometimes show single, internal silt to fine sand laminations that may represent fall overturning.

Figure 2 Time–distance diagram for NAVC varve sections, highlighting sections sampled in this paper. Thin lines show the age range and location of varve sections in central New England that are correlated with the NAVC. Circles and squares highlight basal or ice-proximal varves, respectively, that indicate the location of the ice margin in a given varve year. Thick lines at labeled sections indicate portions of these sections analyzed in this study. As in Fig. 1, the 930PN section in New York is not shown here. The calendar year BP timescale at left uses the value of 20 820 years for the NAVC year–cal yr BP offset from Thompson et al. (2017); see discussion in Sect. 1.2.

1.2 Existing calendar year calibration of the NAVC

Although the NAVC is continuous for its entire length, the lakes in which varves were deposited drained ∼ 12 000 yr BP, so the varve chronology does not have an absolute connection to the calendar year timescale. The majority of the NAVC that lies near or below the glacial–paraglacial transition is replicated in multiple sections (Fig. 2) and is therefore believed to have zero or negligible counting uncertainty; a number of spurious or missing varves at individual sections have been identified and removed from the chronology (Ridge et al., 2012). The uppermost part of the NAVC consists of thin paraglacial varves that only occur at one section (Newbury; see Figs. 1, 2), and Ridge and Toll (1999) estimated that counting uncertainty in the unreplicated part of the Newbury section above ∼ NAVC 7500 is ($+\mathrm{35}/-\mathrm{20}$) years at the top of the section. However, this counting uncertainty is not relevant for the present paper because we will not consider any ice core–NAVC correlations for varves above the glacial–paraglacial boundary. Therefore, if we assume that the NAVC is linear with respect to true calendar years, the calibration of the NAVC to the true calendar year timescale can be represented by a single value for the offset between NAVC years and calendar years. The offset throughout this paper is defined as the age in years BP of NAVC year zero, wherein “years BP” means “years before 1950”, as typically defined for purposes of radiocarbon age calibration. NAVC years are positive and increasing towards the present, and calendar years BP are positive and decreasing towards the present, so the offset is equal to the sum of the NAVC year of any varve and the calendar year BP in which that varve was deposited. Although previous publications (e.g., Ridge et al., 2012; Thompson et al., 2017) variously report estimates for (NAVC–calendar years BP) and (NAVC–calendar years before 2000) offsets, in this paper we use only the NAVC–calendar year BP offset to minimize confusion. Later in this paper, we will also correlate the NAVC with ice core records on the GICC05 timescale (Andersen et al., 2006). Again, throughout this paper we apply the GICC05 timescale in years BP rather than years before 2000 as is sometimes used elsewhere. As the GICC05 timescale is not exactly linear with the true calendar year timescale because of imprecision in ice core layer counting (Adolphi et al., 2018, and references therein), NAVC–cal yr BP and NAVC–GICC05 yr BP offsets may not be exactly the same, and we will use both of these terms as appropriate throughout this paper.

The value of the NAVC–cal yr BP offset has been determined by fitting radiocarbon ages of plant macrofossils associated with individual varves to the INTCAL09 (in Ridge et al., 2012) and subsequently INTCAL13 (in Thompson et al., 2017) calibrated radiocarbon timescales. Thompson et al. (2017), using the INTCAL13 calibration, obtained a best estimate for the NAVC–cal yr BP offset of 20 820 years. The uncertainty in this estimate is not well quantified because residuals between the measured ¹⁴C ages and those predicted by the best-fitting offset to INTCAL13 scatter significantly more than expected from measurement uncertainty and are not normally distributed (see discussion in Ridge et al., 2012, and also in Sect. 4.2.4). Ridge et al. (2012) estimated the uncertainty in the calibrated value of the offset to be on the order of ∼ 200 years and also interpreted the distribution of the residuals to indicate that the most likely source of excess scatter was small amounts of sample contamination by postdepositional carbon. If true, this would imply that the true value of the offset is ∼ 100–200 years higher than the 20 820-year estimate.

Both Ridge et al. (2012) and Thompson et al. (2017) proposed that the estimate of the NAVC–cal yr BP offset could be further refined by aligning prominent climate variations evident in Greenland ice core records with changes in varve thickness and/or ice margin position recorded in the NAVC. For sections of the NAVC near 14 000 yr BP, this alignment yielded a NAVC–GICC05 offset of 20 820 years, which is the same as the best-fitting offset inferred from matching ¹⁴C data to INTCAL13. This result implies a zero offset between INTCAL13 calendar years and GICC05 years at 14 000 yr BP, which is within the range of existing estimates (Adolphi et al., 2018). For sections of the NAVC between 16 000 and 18 000 BP, their alignment of NAVC and ice core climate proxies yielded a NAVC–GICC05 offset of 20 750 years, which may be consistent with the observation that GICC05 likely undercounts ice layers in the 13 000–20 000 yr BP interval (Andersen et al., 2006; Rasmussen et al., 2008; Adolphi et al., 2018). However, as noted above, radiocarbon calibration of the NAVC still suggests that the NAVC–INTCAL13 offset may be ∼ 100–200 years larger than 20 820 years. In addition, it is important to note that matching variations in NAVC and ice core climate proxy records relies on the assumption that the responses of δ¹⁸O variations in Greenland ice cores and varve thickness records in New England to climate changes have similar lag and relative amplitude, which has not been independently established. Also, matches based on these correlations were derived by looking for possible correlations that were close to the NAVC–calendar year offset already estimated from the ¹⁴C data. Both ice core and varve records show periodic variability, so multiple matches could be possible in some time ranges.

1.3 Beryllium-10 in lacustrine sediments

Beryllium-10 fluxes to the Earth's surface in polar regions reconstructed from ¹⁰Be concentrations in ice cores have been widely used as a proxy for solar variability on interannual to centennial timescales (e.g., Beer et al., 1990; Yiou et al., 1997; Finkel and Nishiizumi, 1997; Steig et al., 1998), and ¹⁰Be fluxes to marine sediments have been used to infer paleomagnetic field strength changes on timescales from thousands to millions of years (e.g., Simon et al., 2016, and references therein). By comparison, relatively few studies have attempted to identify solar variability on decadal to centennial timescales using ¹⁰Be measurements from lacustrine sediments. Ljung et al. (2007) found that ¹⁰Be fluxes recorded in lacustrine sediments deposited between 900 and 1450 CE were similar to those recorded in Greenland ice cores for the same time period. More relevant to this study, Mann et al. (2012) and Czymzik et al. (2016, 2018) measured ¹⁰Be concentrations in annually laminated modern and late-glacial sediments, respectively, from European sites and concluded that variability in ¹⁰Be fluxes to the lake was attributable to both direct variations in fallout to the lake due to solar variability and variability in environmental factors affecting transport and scavenging of fallout ¹⁰Be in the lake. These authors suggested several approaches to accounting for these factors in reconstructing the ¹⁰Be fallout flux, mainly focused on multiproxy analysis of lake sediment and the use of multivariate correlation analysis to isolate variability in the ¹⁰Be flux that could not be explained by measured environmental proxies. In this work, we attempt to simplify approaches used in those studies by measuring both ⁹Be and ¹⁰Be in lake sediments and applying the following interpretive framework.

Our basic approach is to (i) use ⁹Be concentrations in sediment as a measure of both the delivery of background Be to the lake and any environmental factors affecting Be scavenging efficiency in the lake and then (ii) identify additional variation in the ¹⁰Be flux that cannot be accounted for by these effects and must therefore represent variability in the atmospheric fallout flux. We quantify this with a simple model for Be fluxes to glacial lake sediments. Because we are working with annually laminated sediments, we can directly measure Be fluxes to the sediment (throughout this paper represented in units of $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) as the product of the Be concentration in sediment (atoms g⁻¹) and the mass accumulation rate ($\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) computed from annual layer thickness (cm yr⁻¹) and density (g cm⁻³).

We assume that measured Be fluxes to lake sediments derive from one source of ⁹Be and two sources of ¹⁰Be. First, sediment delivered to the lake either from glacial discharge or landscape runoff contains “inherited” ⁹Be and ¹⁰Be. Inherited ⁹Be includes adsorbed ⁹Be resulting from subglacial and subaerial mineral weathering. Inherited ¹⁰Be could arise from a number of processes, presumably mainly including subglacial recycling of sediment with adsorbed ¹⁰Be from a past history of surface exposure and interaction between subglacial sediment and fallout ¹⁰Be already present in ice sheet ice. We refer to all Be arising from these sources collectively as “inherited” Be.

Because Be is insoluble and effectively adsorbed to sediment at neutral pH, we expect that sediment delivered to the lake by glacial meltwater discharge will include inherited ⁹Be and ¹⁰Be with a characteristic ¹⁰Be∕⁹Be ratio. The factors controlling the ¹⁰Be∕⁹Be ratio of glacial sediment are presumably complex and include the Be concentration of sediment source material, water–rock interactions in the subglacial environment, the concentration of ¹⁰Be in ice, and the rate and extent of basal and surface melting. However, we assume that whatever these processes, subglacial sediment mixing is effective across a significant area of the ice sheet and results in a characteristic ¹⁰Be∕⁹Be ratio in subglacially discharged sediment that can be assumed to be constant, or at least slowly varying, over the 100- to 1000-year timescales represented by our sample sets.

Second, once glacial sediment with a characteristic ¹⁰Be∕⁹Be ratio is delivered to a proglacial lake, additional ¹⁰Be can be supplied to the lake from atmospheric fallout and adsorbed to the sediment during transport or deposition, but ⁹Be cannot. Thus, Be in lake sediment includes inherited ⁹Be, inherited ¹⁰Be, and fallout ¹⁰Be.

With these assumptions, Be concentrations in lacustrine sediment (which we can measure) can be related to the flux of atmospheric fallout ¹⁰Be (which we seek to reconstruct) by

$$\begin{array}{}\text{(1)}& Q_{\mathrm{10},\mathrm{T}}=Q_{\mathrm{9},\mathrm{T}}{R}_{\mathrm{S}}+Q_{\mathrm{10},\mathrm{A}}f,\end{array}$$

where Q_10,T is the total ¹⁰Be flux ($\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) to the sediment, which is the product of the measured ¹⁰Be concentration (atoms g⁻¹) and the mass accumulation rate ($\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$); Q_9,T is the corresponding total ⁹Be flux ($\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; calculated in the same way as Q_10,T); Q_10,A is the flux of atmospheric fallout ¹⁰Be to the surface ($\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$); R_S (nondimensional) is the characteristic ¹⁰Be∕⁹Be ratio of sediment supplied to the lake; and f (nondimensional) is a focusing factor. The two terms on the right-hand side of the equation represent the two possible sources of ¹⁰Be described above: inherited ¹⁰Be (Q_9,TR_S) and fallout ¹⁰Be (Q_10,Af).

The focusing factor f reflects the fact that ¹⁰Be deposition occurs throughout the lake and perhaps an adjoining watershed on the ice sheet and/or the surrounding ice-free landscape, whereas sediment deposition only occurs in a fraction of the lake basin. Thus, in area-normalized units ($\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), the deposition rate of fallout ¹⁰Be to lake sediment (Q_10,Af) is greater than the fallout rate from the atmosphere (Q_10,A) to the watershed. The focusing factor depends on the geometry of the lake and its watershed; for example, Ljung et al. (2007) estimated f≃20 for a small lake, Belmaker et al. (2008) estimated f≃15 for the much larger Lake Lisan, and the results of Aldahan et al. (1999) implied f≃3 for the still larger Lake Baikal. Absent significant changes in the geometry of the lake and watershed, however, f is expected to be constant or changing slowly, so because we are interested in the short-term variability of the atmospheric fallout flux rather than its absolute value, it is sufficient for our purposes to estimate the quantity Q_10,Af. Therefore, given the framework of Eq. (1) and measured total ⁹Be and ¹⁰Be fluxes to the sediment, the primary obstacle to reconstructing atmospheric fallout fluxes of ¹⁰Be is the need for an estimate of R_S, which we discuss in detail later from the perspective of analyzing each of our data sets.

The purpose of the model described by Eq. 1 is to separate environmental effects influencing the delivery of Be to the lake from variations in ¹⁰Be fallout to the lake and watershed. If the assumptions of the model and the estimate of R_S used in applying it are correct, variations in Q_10,Af reconstructed by this method should correctly represent variations in fallout and not in Be transport processes. However, the model is very simple and several complications may be relevant for our data. One is that R_S and f will likely change over time due to (i) retreat of the ice margin that results in an increase in the proportion of sediment supplied by the deglaciated portion of the watershed, (ii) changes in ice-marginal drainage, and eventually (iii) reduction in glacial meltwater and sediment supply in the transition from a glacial to a paraglacial lake. As noted above, however, it is likely that these changes will either be steady, gradual changes associated with ice margin retreat from the lake basin or large instantaneous changes associated with lake inlet or outlet switching or reconfiguration of ice-marginal drainage. It appears unlikely that these processes would be manifested as periodic decadal- to centennial-scale variations that might mimic the solar variability we seek to identify. A second complication is the possible role of watershed processes; if fallout ¹⁰Be were sourced from a significant ice-free watershed surrounding the lake, variations in ¹⁰Be fallout rates might be buffered by water–soil interactions in the watershed and reduce observed fallout variability in the lake (e.g., Czymzik et al., 2016). Finally, the potential importance of ¹⁰Be already present in ice sheet ice that could be released by surface melt is unknown. The expected contribution from typical ¹⁰Be concentrations in ice sheets (order 10⁴ atoms g⁻¹) at moderate melt rates averaged over the ablation zone (order 10 cm yr⁻¹) would be small (order 10⁵ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) compared to the atmospheric fallout flux (order ∼ 1.5×10⁶ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; see Graly et al., 2011). On the other hand, if melt rates averaged over the contributing ablation area reached meters per year, which could be possible during periods of rapid retreat, this contribution might be significant.

2 Methods

2.1 Sample collection

We collected several sample sets from a total of seven varve sections, with various sampling strategies designed to address different aspects of the study (Table 1 and discussion below). Samples were extracted from (i) short cores collected from surface outcrops by cleaning outcrops and driving in short (typically 0.6 m) sections of PVC pipe with a steel cap and sledge hammer as well as from (ii) deep (up to 50 m), continuously sampled hollow-stem auger cores collected in 2007–2009. Cores are archived at Tufts University. All cores were split after collection and partially dried to improve visual identification of layers and facilitate sampling of individual varves. Only cores with minimal coring deformation were sampled. We sampled winter and summer layers of single varves, complete single varves, and continuous sets of contiguous varves by scraping core surfaces and cutting interior sections of the cores to remove any possible surface contamination, then separating varves along partings between layers and scraping top and bottom surfaces to isolate the varve(s) of interest. Each sample was a column or wedge cut such that the cross-sectional area was constant for all varves incorporated in the sample. As in previous work relating to the NAVC (Ridge et al., 2012), we consider a single varve to extend from the bottom of the coarse-grained summer layer to the top of the fine-grained winter layer.

2.1.1 Short biennial records

At two sites, one with glacial varves (the Kelsey–Ferguson, or KF, section; Figs. 1, 2; Table 1) and one with paraglacial varves (the North Haverhill, or NHV, section; Figs. 1, 2; Table 1), we collected continuous sets of samples spanning 80-year periods at 2-year resolution. The purpose of these sample sets is to determine whether decadal-scale solar variability, specifically the ∼ 11-year Schwabe cycle, can be identified in ¹⁰Be fluxes to NAVC sediments (see Sect. 4.1 below).

2.1.2 Long decadal records

At two sites, we collected sets of samples designed to produce records of ¹⁰Be deposition over 1000–1500-year periods with decadal (ca. 25-year) resolution that could potentially be correlated with ¹⁰Be records from ice cores (see Sect. 4.2 below). One set of samples is from glacial varves at the Kelsey–Ferguson (KF) and nearby Glastonbury (GL) sites in central Connecticut (Figs. 1, 2; Table 1) that overlap and together span NAVC years 3658–4675. The second set is a sequence of glacial grading into paraglacial varves from Newbury, Vermont (NEWB; Figs. 1, 2; Table 1), that spans NAVC years 6631–8193.

2.1.3 Winter–summer pairs

To investigate the seasonal distribution of Be deposition, we separately analyzed winter and summer layers from individual varves at four sites (Figs. 1, 2; Table 1): (i) glacial varves from the Kelsey–Ferguson (KF) site; (ii) paraglacial varves from North Haverhill (NHV); (iii) glacial varves at a site in North Hatfield (NHT), MA; and (iv) glacial varves from a section in the Mohawk Valley of central New York (930PN) that is not correlated with the NAVC but was deposited approximately 17 000 yr BP. Varves at 930PN are unusual among NAVC sites in that their parent material is limestone, so they are mostly carbonate.

2.1.4 Ice-proximal sediment and underlying till

To investigate inherited Be concentrations in subglacial sediment that serves as the parent material for varved lake sediments, we analyzed basal ice-proximal varves and underlying diamicton at South Hadley (SHD) and North Hatfield (NHT), MA, from hollow-stem auger cores that penetrated the base of these varve sections. The basal varves are relatively coarse, 2–6 cm-thick varves deposited at the bottom of these sections, presumably close to the ice margin. The diamicton at SHD is a stony subglacial till, whereas that at NHT mixes till and glaciolacustrine sediment and likely represents a subaqueous debris flow near the ice sheet grounding line.

2.2 Grain size measurements

We measured grain size distributions for all samples using a combination of mechanical sieving (for grain sizes larger than 4 Φ) and pipette analysis for smaller grain sizes (e.g., Lewis and McConchie, 2012). Throughout this paper we represent grain size using the logarithmic phi scale of Krumbein (1934), which is commonly used in sedimentology: Φ is the negative base-2 logarithm of the grain diameter in millimeters, so larger values of Φ represent smaller grain sizes. We halted the pipette analysis at 11 Φ (∼ 0.5 µm), but significant fractions (up to 60 %) of all samples were finer than this (Fig. 3). Thus, to calculate mean grain size and sorting parameters (see, e.g., Boggs, 2014) that require the complete distribution, we fit a lognormal distribution to the observed data for each sample and used the fitted lognormal distribution to compute mean grain size and sorting (Fig. 3).

Figure 3 Grain size and sediment density measurements. (a) Measured cumulative grain size distributions for summer layers (red), winter layers (blue), and complete single or multiple varves (gray). Pipette analysis stopped at 11 Φ. (b) Example of extrapolating measured grain size data for a representative sample (NAVC 7068-69 from the NHV section) using a fitted lognormal distribution. (c) Relationship between sediment dry density and grain size. The linear approximation shown by the black line (which is fit to all data excluding three low outliers that likely result from failure to completely fill the reference volume with stiff glacial sediment) is used in subsequent flux calculations.

2.3 Sediment density

Measurements of sediment density are required to interpret a Be concentration by weight (atoms g⁻¹) as a Be flux to the sediment ($\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). We measured the dry density of representative samples of glacial varves at the Kelsey–Ferguson (KF) site, paraglacial varves from the North Haverhill (NHV) site, and separate summer and winter layers from North Hatfield (NHT) by pressing a cylinder of known volume into wet sediment and drying and weighing the resulting sample (Fig. 3; Table S2 in the Supplement). This procedure has several potential sources of inaccuracy, primarily due to the stiff and friable nature of the sediment, which causes difficulty in completely filling the standard volume without deforming or compressing the sediment or introducing voids. It is important to note that the purpose of this study is to identify variations in Be flux to the sediment, not to make an accurate estimate of the absolute flux, so from this perspective inaccuracies in density measurements are not necessarily significant. On the other hand, if we assumed a constant density when, in fact, periodic variations in density were present, we could infer spurious Be flux variations. Density variations, in turn, are most likely to be the result of variations in grain size and therefore in porosity and the relative proportion of silt and clay, and our measurements show this effect (Fig. 3). Thus, for subsequent Be flux calculations we compute the density of all samples from mean grain size measurements using the linear relationship shown in Fig. 3, with an uncertainty derived from the standard deviation of the residuals (±0.05 g cm⁻³).

2.4 Beryllium-9 by ICP-OES

We measured concentrations of adsorbed ⁹Be using the procedure described in Greene (2016), which involves drying and powdering sediment samples, leaching in 6M HCl for 24 h, then diluting to a suitable concentration for analysis and measuring the Be concentration with a JY Horiba Optima inductively coupled plasma optical emission spectrometer (ICP-OES). Greene (2016) discusses step-leaching experiments and other data used for quality control by assessing the reliability of the procedure. This procedure extracts only adsorbed Be and not structural Be contained in Be-bearing minerals (e.g., beryl). Although the nominal internal precision of the ICP-OES measurement inferred from replicate analysis of liquid standards was typically 1 %–2 %, replicate analyses of splits of the powdered sediment samples showed a total uncertainty of 3.5 %. Thus, we assume that all ⁹Be measurements have 3.5 % uncertainty.

2.5 Beryllium-10 by accelerator mass spectrometry

We measured ¹⁰Be concentrations in sediment samples at the University of Vermont (UVM) using an adaptation of the potassium bifluoride fusion method of Stone (1998). We dried sediment samples, homogenized them by powdering in a Spex shatterbox mill, and subsampled 0.5 g of the homogenized sediment. We mixed the sample with KHF₂ and Na₂SO₄ and spiked the mixture with 0.4 mg of one of several Be carriers prepared from deep-mined beryl and having ¹⁰Be∕⁹Be in the range 2–$\mathrm{8}\times {\mathrm{10}}^{-\mathrm{16}}$. We then extracted Be by fusion of the sample, leaching of the fused sample in water, precipitation of excess K as KClO4, and precipitation of Be hydroxide. We modified the published method by adding an ion exchange step in which we redissolved Be hydroxide in dilute HCl, loaded it onto cation exchange resin (Bio-Rad AG50W X8), eluted B and other cations in dilute HCl, recovered Be in 6M HCl, evaporated to dryness, redissolved in dilute HCl, and precipitated Be hydroxide. Finally, we converted Be hydroxide to BeO for Be isotope ratio measurement by accelerator mass spectrometry (AMS) at the Scottish Universities Environmental Research Centre (SUERC) in East Kilbride, Scotland (Xu et al., 2015). For these analyses, we sputtered samples and isotope ratio standards in the ion source for the same amount of time in order to suppress any effect of variations in emittance during sputtering. All ¹⁰Be measurements are normalized to the NIST SRM4325 standard material with an assumed ¹⁰Be∕⁹Be ratio of $\mathrm{2.79}\times {\mathrm{10}}^{-\mathrm{11}}$, which is equivalent to the 07KNSTD standardization of Nishiizumi et al. (2007). We processed samples in batches of 16, each containing one process blank and at least one replicate of a sample that was also analyzed in a different batch (see below). The process blank for Be fusion extraction at UVM utilizes a sample of fine-grained fluvial sediment chosen because of its naturally low bulk ¹⁰Be concentration and subsequently powdered, homogenized, and etched in dilute HNO₃ in a heated sonic bath for an extended period to remove adsorbed Be. We processed 0.5 g aliquots of this material in the same way as the samples. Complete carrier and process blanks contained between 120 000 and 310 000 atoms ¹⁰Be, representing 0.1 %–0.7 % of total atoms measured in samples.

In contrast to the leaching procedure used to measure ⁹Be concentrations, the total fusion method extracts both adsorbed ¹⁰Be and mineral-bound ¹⁰Be that could result from in situ production of cosmogenic ¹⁰Be in any fraction of the sediment that previously experienced surface exposure. However, concentrations of in-situ-produced cosmogenic ¹⁰Be are expected to be negligible for the purposes of this study because in situ cosmogenic production rates at moderate elevations (order 10 $\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) are insignificant compared to deposition rates of atmospherically produced ¹⁰Be (order 10⁶ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and because NAVC sediments are mostly subglacially derived and therefore not previously exposed to the surface cosmic-ray flux. In-situ-produced cosmogenic ¹⁰Be concentrations in glacial sediments associated with the Laurentide Ice Sheet, where measured (e.g., Balco et al., 2005), are on the order of 10 000 atoms g⁻¹, which is negligible compared to total ¹⁰Be concentrations of order 10⁸ atoms g⁻¹ measured in this study. We therefore assume that the technical inconsistency between our measurements of leachable ⁹Be and total ¹⁰Be is insignificant and disregard it.

To quantify total measurement uncertainty and as a general quality control measure, we prepared large quantities of two representative samples to be analyzed repeatedly throughout the project. The first of these consisted of varves NAVC 3569-3570 from the Kelsey–Ferguson site, collected as part of the 80-year biennial record described above, which we refer to here by the internal UVM sample name BD3. When the BD3 sample was exhausted after 15 measurements, we prepared a second sample for replicate analysis (BD-STAN) by mixing and homogenizing aliquots of several other glaciolacustrine sediment samples from this study. The BD-STAN sample was subsequently analyzed six times.

The purpose of chemical purification of Be for AMS analysis is to produce a sample of BeO that is sufficiently pure to produce a Be ion beam of the same intensity as produced by pure reagent BeO used as an isotope ratio standard. Poor or inconsistent Be yield from chemical purification, or failure to consistently remove contaminants from the BeO target, can result in Be ion beam currents that are low and/or inconsistent compared to those from standards (see Hunt et al., 2008). This is undesirable because (i) low beam currents reduce measurement precision by reducing ¹⁰Be count rates, and (ii) inconsistent beam currents can introduce systematic errors in isotope ratio measurements due to current-dependent variations in beam emittance, ion beam transport, or detector efficiency. Although the goal of AMS beamline and detector setup and adjustment is to remove or minimize any beam current dependence of ratio measurements, it is difficult to verify whether this has been achieved for an arbitrarily large range in beam current, so instrument tuning is not a complete substitute for consistent sample preparation.

AMS analysis of Be samples prepared in the first eight chemistry batches showed low and inconsistent beam currents (Fig. 4; mean and standard deviation of sample beam currents relative to standards for these batches were 0.63 ± 0.25). We found that this effect was related to trace Cl remaining in samples after hydroxide precipitation, so we further modified the Be extraction procedure by redissolving Be hydroxide in dilute HNO₃ and reprecipitating. Be beam currents from samples in 17 chemistry batches after this modification were higher and less variable (Fig. 4; mean and standard deviation of relative sample beam currents after the modification were 0.87 ± 0.16).

Figure 4 Be beam currents during AMS analysis at SUERC relative to Be isotope ratio standards. Each symbol shows the average beam current during analysis of a single sample, normalized to the average current of all primary standards analyzed during the same measurement period. Red lines and pink bars show means and standard deviations for samples processed in each chemistry batch. UVM batch ID numbers are assigned consecutively, so batches are shown in the order they were processed. Only samples of NAVC sediments that are part of this study are included in this plot; some batches included unrelated samples, which are not shown.

Although this change to the sample preparation method significantly improved beam current performance, the complete data set displays high variability in Be currents. Thus, we used the results of replicate analyses, which spanned a range of Be currents, to look for and quantify any beam current dependence in the AMS results. In both replicate analyses and in analyses of groups of samples with similar ¹⁰Be concentrations, we observed a relationship between Be beam current and apparent ¹⁰Be concentrations derived from AMS-measured ¹⁰Be∕⁹Be ratios (Fig. 5). Although apparent concentrations do not display a significant correlation with beam current for samples with beam currents similar to standards (between ca. 80 % and 110 % of standard currents), our data include a much wider range of beam currents, and when all data are considered, significant correlations are present (Fig. 5). We observed this effect to be persistent through many AMS measurement periods during a 3-year period and similar to the current dependence observed for many replicate measurements of in-situ-produced ¹⁰Be in the CRONUS-N quartz standard measured at the SUERC accelerator by Bierman et al. (2017) as well as for analysis of the CoQtz-N quartz standard at an anonymous AMS facility described by Binnie et al. (2019, their Fig. 3).

Figure 5 Observed relation between ⁹Be beam current relative to standards and the measured ¹⁰Be∕⁹Be ratio, here expressed as the nuclide concentration instead of the ratio to account for small differences (ca. 5 %) in sample and Be carrier masses between replicate analyses (see text). (a) 15 replicate analyses of the BD3 sample (see text). (b) Replicate analyses of samples from 80-year biennial records at the KF and NHV sites (see text above); each of these samples was measured at least twice, and light-colored lines connect replicate analyses. Note that true ¹⁰Be concentrations vary among samples in this data set, but a significant correlation with beam current is present regardless, and nearly all replicates of individual samples display a positive slope in this diagram. In the left-hand two panels, symbol color-coding indicates the AMS measurement session in which each analysis took place. Sessions 1–2 were measured in 2012, 3–5 in 2013, and 6–9 in 2014. (c) Equivalent data for in-situ-produced ¹⁰Be in the CRONUS-N quartz standard (Jull et al., 2015) prepared in the UVM chemistry lab and measured on the SUERC accelerator between 2013 and 2017, as reported by Bierman et al. (2017).

We conclude from these observations that apparent measured ¹⁰Be concentrations for our samples that had anomalously low or high AMS beam currents are inaccurate due to beam current dependence of the measured isotope ratios. Lacking a physical model for the observed beam current dependence, we devised the following procedure to correct for this bias. If (C∕C_S) is the beam current for a sample relative to the average current for standards measured simultaneously, R_M is the apparent ¹⁰Be∕⁹Be ratio observed for the sample, and R_C is a nominally correct ¹⁰Be∕⁹Be ratio that would be observed for $(C/C_{\mathrm{S}})=\mathrm{0.9}$, then we can define an empirical linear correction factor S such that $R_{\mathrm{M}}/R_{\mathrm{C}}=\mathrm{1}-S(\mathrm{0.9}-C/C_{\mathrm{S}})$. In the case in which blank corrections are small to negligible (as is the case for these data; see above) the ratio of apparent and corrected nuclide concentrations N_M∕N_C, where N_M and N_C are nuclide concentrations in atoms g⁻¹, is equal to R_M∕R_C. Thus, to account for small variations in sample and ⁹Be carrier masses among replicate samples, we correct apparent concentrations instead of ratios using $N_{\mathrm{M}}/N_{\mathrm{C}}=\mathrm{1}-S(\mathrm{0.9}-C/C_{\mathrm{S}})$. 0.9 is chosen as the normalizing value because it is the modal value of (C∕C_S) for our complete data set, so it minimizes required corrections. Choosing a normalizing value not equal to 1 might introduce a small systematic bias in measured ¹⁰Be concentrations relative to the primary measurement standards, but that would not affect any of the results in this study because we are concerned with variability in the ¹⁰Be flux and not with a precise measurement of the absolute flux.

Figure 6 Correction for current dependence bias. Uncorrected data shown by open circles are the same data shown in Fig. 5. Corrected concentrations shown by filled circles are connected by light lines to corresponding uncorrected concentrations; error bars are omitted for clarity. For measurements with beam currents comparable to standards, corrections are minimal and similar in size to estimated measurement uncertainties; corrections are most significant for measurements with extreme beam current values. The correlation coefficients and probabilities shown at upper left in all panels are for the corrected concentrations.

We estimate S by choosing the value that minimizes scatter among replicate analyses of the same sample. The sample with the largest number of replicates (BD3; see discussion above) yields S=0.164. The CRONUS-N data of Bierman et al. (2017) yield S=0.185, which highlights the apparent consistency of this effect across multiple samples, chemical preparation methods, and AMS measurement periods. Figure 6 shows the effect of applying this correction procedure with S=0.174 (the average of the above two values) for the data in Fig. 5. Corrected data are uncorrelated with relative beam currents. Note that this correction scheme has the property that corrections for measurements with beam currents comparable to standards (C∕C_S between ca. 0.8 and 1.1, within which range no significant correlation between the current and ratio is observed) are minimal and comparable to individual measurement uncertainties; corrections only become significant for samples with extreme relative currents. Figure 7 further highlights this point for the 80-year biennial record from North Haverhill (NHV), which was entirely measured in duplicate. The correction procedure improves agreement between the two sets of measurements primarily by correcting a minority of data that were biased low due to anomalously low beam currents without having a significant effect on the majority of the measurements. For the remainder of this paper, we correct all ¹⁰Be concentrations for beam current bias using S=0.174. Although presumably this estimate of S has some uncertainty, a formal estimate of the uncertainty is not necessary because, as discussed below, we derive total measurement uncertainties for ¹⁰Be concentrations from the statistics of replicate analyses after correction for beam current bias rather than by propagation of errors through the correction process.

Figure 7 Effect of current dependence bias correction on agreement between replicate analyses of NHV biennial samples. All data are assigned 4 % uncertainties as discussed below.

Nominal measurement uncertainties on corrected ¹⁰Be concentrations include: (i) AMS measurement uncertainties derived from machine stability and counting statistics (∼ 1 %) and replicate analyses of secondary isotope ratio standards during AMS measurement periods (∼ 1 %); (ii) uncertainty in ⁹Be carrier concentrations (<1 %); (iii) blank uncertainty (<0.2 %); and (iv) uncertainty in S as discussed above (unknown). Combined, these indicate measurement uncertainties no less than 1.5 %–2 % for most samples. Scatter in measured ¹⁰Be concentrations in replicate aliquots of homogenized sediment samples, however, is significantly higher, which indicates that propagation of known uncertainties would underestimate the true measurement uncertainty. Figure 8 shows the reproducibility of replicate analyses of a total of 84 samples. The mean relative standard deviation of all sets of replicate analyses is 4 %, that for BD3 is 3.5 % (for 15 measurements), and that for BD-STAN is 4.1 % (for 6 measurements). We conclude that the true uncertainty in a typical measurement from our sample set is near 4 %, which likely reflects the known measurement uncertainties enumerated above, errors introduced by the correction for current dependence, and imperfect homogenization of sediment samples during preparation. Thus, we compute total uncertainties for all ¹⁰Be concentration measurements as follows. First, we assume that each individual measurement has 4 % precision. Second, to incorporate the principle that repeated measurements of the same sample should improve total measurement precision, for samples analyzed more than once, we calculate a weighted mean and standard error (although, as all individual measurements are assigned 4 % uncertainty, they are weighted equally). Summary concentrations and uncertainties calculated via this approach appear in the Supplement and are used henceforth throughout the paper.

Figure 8 Relative standard deviation for analyses of replicate splits of 84 sediment samples that were analyzed more than once. Data have been corrected for current dependence bias.

3 Results and discussion I: general framework and systematics

Beryllium-10 concentrations for all NAVC varves (Fig. 9;Table S1) are in the range 0.5–1.5×10⁸ atoms g⁻¹, which is typical of concentrations observed in lacustrine sediments elsewhere (typically 1–10×10⁸ atoms g⁻¹; e.g., Ljung et al., 2007; Czymzik et al., 2016). Total leachable Be concentrations in sediments are 0.6–1.7 ppm, so ¹⁰Be∕⁹Be ratios are in the range 1–$\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$. As expected from the inverse relationship between grain size and surface area (e.g., Pavich et al., 1984), the dominant control on adsorbed Be concentrations is grain size: both ⁹Be and ¹⁰Be concentrations are inversely correlated with mean grain size at all sites (Fig. 10). However, ¹⁰Be∕⁹Be ratios are not correlated with grain size, which is important because it supports the hypothesis used in formulating Eq. (1) that inherited Be adsorbed to sediment has a characteristic ¹⁰Be∕⁹Be ratio.

Figure 9 Mean grain size and ¹⁰Be–⁹Be concentrations for all samples, shown in stratigraphic order by NAVC year. Seasonal samples from site 930PN, which is not correlated with the NAVC, are not shown. The boundary between glacial and paraglacial varves is between NAVC 7020 and 7200 but is transitional at many sites and also not isochronous between sections. Error bars are not shown for clarity.

Variations in mean grain size, and thus in Be concentrations, at our sample sites reflect a balance of two primary effects. Initially, retreat of the ice margin from a site and therefore an increase in distance from the glacial sediment source will cause an upsection decrease in mean grain size. Subsequently, water depth decreases at the site due to sediment infilling of the lake and glacioisostatic rebound relative to an ice-distal outlet. Shallowing increases the energy of the depositional environment because of increased exposure to surface waves and because the decreasing cross-sectional area of the lake requires that current velocities increase to maintain constant discharge, and it also leads to an upsection increase in grain size. An additional possible effect that may be relevant at some sites is that as the ice margin retreats, sediment is increasingly sourced from the surrounding deglaciated landscape rather than direct ice sheet discharge, and paraglacial runoff sediment is likely coarser than subglacial sediment.

Figure 10 Correlation analysis for grain size and Be isotope concentrations and ratio. The correlation coefficient and p value are shown in red when data are correlated at 95 % or better confidence. ⁹Be and ¹⁰Be concentrations are always significantly correlated with mean grain size. However, the ¹⁰Be∕⁹Be ratio is not. Note that Φ size is inverse to dimensional grain size, so a positive correlation between Be concentration and Φ is a negative correlation between Be concentration and grain size.

The long record at Newbury that we sampled with decadal resolution (NAVC 6631–8193; Fig. 9) shows these effects. Grain size decreases and Be concentrations correspondingly increase in the initial part of the record (NAVC 6600–6900) as the ice margin recedes from the site, but subsequently lake shallowing becomes more important and these trends reverse (NAVC 6900–8100). The long records at the Kelsey–Ferguson and Glastonbury sites (NAVC 3658–4669; Fig. 9), on the other hand, begin in relatively ice-distal varves because these sites were deglaciated at about NAVC 2800. Thus, grain size increases and Be concentrations decrease throughout, following the pattern in the upper half of the Newbury section.

Measurements of Be concentrations in separate summer layer – winter layer pairs from individual varves (Figs. 9, 10) highlight the grain size dependence of Be concentrations: finer-grained winter layers have much higher Be concentrations than summer layers, but ¹⁰Be∕⁹Be ratios are not significantly different. In addition, summer and winter data display the same grain size–concentration relationship as full-year data (Fig. 10). Although the winter and summer samples define the end members of the distribution in grain size–concentration space, seasonal and full-year data lie on the same trend.

Samples of ice-proximal glaciolacustrine sediment and underlying diamict (see Sect. 2.1.4) have similar Be concentrations and isotope ratios as more ice-distal glaciolacustrine sediment in similar stratigraphic positions (compare the green symbols in Fig. 9 to other data), with one exception. The exception is a subglacial till (from the SHD section; green square in Fig. 9). In contrast to all other samples, (i) the properties and stratigraphic position of this sample indicate that it has no glaciolacustrine input and has most likely not interacted with lake water, and (ii) this sample has much lower ⁹Be and ¹⁰Be concentrations. Lower Be concentrations are consistent with the relatively large mean grain size (1.1 Φ compared to 8–11 Φ for glaciolacustrine samples, although with much poorer sorting). However, the much lower ¹⁰Be∕⁹Be ratio ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$) in this sample than in glaciolacustrine sediments (∼ 1–$\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$; see Fig. 9) is not necessarily expected and may indicate that the majority of inherited ¹⁰Be in glaciolacustrine sediments is not reflective of recycling of previously exposed sediment, but it instead may derive from interaction within the subglacial water system between the till and ¹⁰Be that is present in ice and released by melting.

4 Results and discussion II: biennial and decadal data series

4.1 Biennial series

The purpose of sampling two 80-year varve sequences at 2-year resolution is to determine whether or not the 11-year Schwabe solar cycle is present in ¹⁰Be flux records. This cycle is evident in ice cores (e.g., Steig et al., 1998) and some lake sediments (e.g., Mann et al., 2012). Recognizing this cycle in the ¹⁰Be flux to NAVC sediments would be strong evidence that the NAVC does, in fact, contain a record of global variations in ¹⁰Be production and fallout.

4.1.1 Characteristics of biennial series

An 80-year segment of the paraglacial varve sequence at North Haverhill (NHV; Fig. 11) has a mean varve thickness of 0.3 cm, corresponding to mass accumulation rates in the range 0.25–0.75 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. As discussed above, ¹⁰Be and ⁹Be concentrations are correlated with mean grain size (Figs. 10, 11). Multitaper spectral analysis (Fig. 11; Dettinger et al., 1995; Ghil et al., 2002) does not show significant spectral power in the 11-year band for any of the observed data series except mean grain size, for which a 12-year peak is significant at 90 % confidence.

Figure 11 Sedimentological data and Be concentrations for an 80-year record with 2-year resolution from paraglacial varves at the North Haverhill (NHV) site. The mass accumulation rate (MAR) is calculated from varve thickness and the density–grain size relationship shown in Fig. 3. Plots at right show results of multitaper method (MTM) spectral analysis of the observational data series using the SSA-MTM Toolkit (Dettinger et al., 1995; Ghil et al., 2002), highlighting frequencies in the vicinity of the 11-year Schwabe solar cycle. Blue and red lines show the 50 % and 90 % significance thresholds, respectively, relative to an AR(1) red noise estimate, and the shaded area highlights the 11-year solar cycle band.

The second 80-year sequence is from glacial varves at the Kelsey–Ferguson (KF) site (Fig. 12) and has a mean varve thickness of 0.85 cm and mass accumulation rates of 0.5–2.5 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, more than twice that for the paraglacial varves at NHV. Although all the data from KF (shown in Fig. 10) display significant correlations between grain size and both ⁹Be and ¹⁰Be concentrations, in this subset of the data, ⁹Be concentrations are significantly correlated with mean grain size (r=0.87; p<0.01) but ¹⁰Be concentrations are not (r=0.18; p=0.28). None of the observational data series show significant spectral power in the 11-year band.

Figure 12 Sedimentological data and Be concentrations for 80-year record with 2-year resolution from glacial varves at Kelsey–Ferguson (KF) site, with MTM spectra. Plot elements are the same as described in Fig. 11. For the purposes of spectral analysis, missing grain size and ¹⁰Be data for two samples were filled by linear interpolation between adjacent measurements.

4.1.2 ¹⁰Be systematics in biennial series

Be isotope fluxes to the sediment can be derived from the data in Figs. 11 and 12 and are shown in Figs. 13, 14, and 15 for the KF and NHV biennial data. As discussed above, the Be flux (in $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) is computed by multiplying Be isotope concentrations (atoms g⁻¹) and mass accumulation rates ($\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Following Eq. (1), the ⁹Be flux is inherited ⁹Be adsorbed to sediment entering the lake, whereas the total ¹⁰Be flux includes both inherited and fallout components.

We are interested in reconstructing the fallout ¹⁰Be flux (Q_10,Af in Eq. 1), and here we attempt to do this by applying Eq. (1) with the following additional assumptions. First, we assume that the inherited ¹⁰Be∕⁹Be ratio R_S is constant throughout the time series, which, as discussed above, appears likely over short periods of time. Second, we assume that the fallout flux Q_10,Af is variable but approximately symmetrically distributed around a mean value. Although atmospheric production and transport processes are nonlinear with respect to solar variability, assuming a symmetric distribution is likely adequate for short-term variability in the fallout flux that is periodic and has relatively low magnitude. For completeness, we also assume that there is no common forcing of fallout variations and environmental factors affecting ⁹Be delivery to the lake, which is supported by the fact that Rittenour et al. (2000) did not observe 11-year-band spectral power in NAVC varve thickness records.

Figure 13 Plots showing regression of ¹⁰Be against ⁹Be fluxes for NHV and KF biennial data. Following Eq. (1), the y intercept should give the mean ¹⁰Be fallout flux. It is evident that, although this is above zero (as it should be), it is small compared to total ¹⁰Be fluxes and comparable in magnitude to measurement uncertainty and scatter around the line. The dotted lines show a 68 % confidence interval for the regression derived from a Monte Carlo simulation.

With these assumptions, Eq. (1) predicts that linear regression of (Q_9,T, Q_10,T) pairs will yield a line having a slope equal to R_S and y intercept equal to the mean fallout flux. Sample-to-sample variation in the fallout flux is expressed as scatter around the line. Q_9,T and Q_10,T are, in fact, linearly related for both biennial records (Fig. 13). Linear regression for NHV yields $R_{\mathrm{S}}=\mathrm{1.539}\pm \mathrm{0.077}\times {\mathrm{10}}^{-\mathrm{9}}$ and mean(Q_10,Af) $=\mathrm{3.0}\pm \mathrm{1.9}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The regression for KF yields $R_{\mathrm{S}}=\mathrm{1.252}\pm \mathrm{0.054}\times {\mathrm{10}}^{-\mathrm{9}}$ and mean(Q_10,Af) $=\mathrm{3.0}\pm \mathrm{4.9}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. However, these results imply several complications. First, although intercepts are positive for both data sets, neither are distinguishable from zero at high confidence, and both imply that the fallout ¹⁰Be flux is a small fraction (< 10 %) of the total ¹⁰Be flux at these sites. This is not favorable from the perspective of reconstructing solar-cycle variability in ¹⁰Be fallout fluxes because the expected magnitude of variability in ¹⁰Be fallout associated with the 11-year cycle is ∼ 10 %–20 % (Lal and Peters, 1967; Steig et al., 1998). A 10 %–20 % variation in 10 % of the total implies variability in total measured ¹⁰Be fluxes no greater than 1 %–2 %, which is smaller than our measurement uncertainties of 3 %–4 %. Second, in agreement with this prediction, Fig. 13 shows that measured ¹⁰Be fluxes at both sites are indistinguishable from the regression line at the level of measurement uncertainties. In other words, from this analysis alone we cannot reject the hypothesis that there is zero variation in the ¹⁰Be fallout flux and deviations from the linear relationship are entirely due to measurement errors. Both of these observations indicate that relatively high concentrations of inherited ¹⁰Be and relatively low concentrations of fallout ¹⁰Be result in a poor signal-to-noise ratio for 11-year-period variations in fallout ¹⁰Be and a low likelihood of detecting these variations.

Figure 14 Be flux estimates for an 80-year biennial record at NHV. Total fluxes are calculated using the grain size–sediment density relationship shown in Fig. 3. ¹⁰Be fallout fluxes are estimated from Eq. (1) using two different values of R_S: the value inferred from the linear regression in Fig. 13 and a second value chosen to yield higher apparent fallout fluxes and highlight the fact that apparent time variations in the flux are not strongly dependent on the choice of R_S. Error bars reflect the propagation of uncertainty on Be isotope concentrations but not on mass accumulation rates, so it may be underestimated. Plot elements for MTM spectra are the same as in Fig. 11.

Figures 14 and 15 show estimates of the ¹⁰Be fallout flux for the biennial data sets using the estimates of R_S derived from the linear regressions in Fig. 13. Note that because the scatter of data around the regression lines has a similar magnitude as the mean fallout fluxes estimated from the y intercept, this calculation results in some values for the fallout flux that are less than or indistinguishable from zero. In addition, this calculation yields fallout fluxes that scatter by ∼ 60 % around the mean, which is substantially more than expected for the 11-year solar cycle. These observations suggest that the regression method may be underestimating the ¹⁰Be fallout flux. A systematic underestimate might not matter if we are only concerned with whether or not periodic variability exists, but to account for this possibility, we repeated the calculation with values for R_S that are lower than inferred from the regressions in Fig. 13 and chosen to yield mean fallout flux estimates that all exceed zero (Figs. 14, 15). Although there is no strong basis for choosing these values, they yield estimates of the ¹⁰Be fallout flux that have less relative variability and highlight the property of Eq. (1) that although the choice of R_S affects the mean estimated ¹⁰Be fallout flux and therefore the relative magnitude of the variability around the mean, it has a minimal effect on its pattern of variability. This, in turn, is important because it implies that an inaccurate estimate of R_S is most likely not an obstacle to spectral analysis of short-term variability or correlation with ice core records. On the other hand, the relative magnitude of variability in the reconstructed fallout fluxes is strongly dependent on the assumed value of R_S. As our estimate of R_S depends on several assumptions and may be unreliable, the magnitude of reconstructed fallout variability is also likely unreliable. This, in turn, is potentially important in evaluating whether or not reconstructed ¹⁰Be fallout variations from NAVC sediments are consistent with independent records or with physical limits on the possible variability in atmospheric ¹⁰Be production and fallout, and we discuss it in more detail later in Sect. 4.2. From the perspective of the biennial records at NHV and KF, however, this issue is not likely to be relevant because, as highlighted in Figs. 14 and 15, propagated measurement uncertainties in the estimates of the ¹⁰Be fallout flux, especially for the paraglacial NHV data, have the same scale as the short-term variability in the estimates. This supports the argument that the unexpectedly large magnitude of short-term variation in the ¹⁰Be fallout flux is not the result of incorrectly estimating R_S, but indicates that variability in the true fallout flux is smaller than measurement uncertainty in the reconstructed fallout flux and that the majority of apparent variability derives from measurement errors.

Figure 15 Be flux estimates for the 80-year biennial record at KF. Total fluxes are calculated using the grain size–sediment density relationship shown in Fig. 3. Fallout fluxes are estimated from two different assumed values of R_S; see text for details. Error bars reflect propagation of uncertainty on Be isotope concentrations but not on mass accumulation rates. Plot elements for MTM spectra are the same as in Fig. 11.

Spectral analysis of reconstructed ¹⁰Be fallout fluxes for both biennial records does not display significant spectral power in the 11-year band (Figs. 14, 15). However, results of spectral analysis highlight that the choice of R_S has limited importance for spectral analysis of the records; although the significance level of individual peaks varies somewhat with R_S, power spectra are similar for different values. For example, R_S derived from the regression analysis at NHV (Fig. 14) yields a 14-year peak that is significant at 90 % confidence; a lower value of R_S yields the same peak but at lower significance. At KF (Fig. 15), on the other hand, a lower value of R_S slightly increases the significance level of peaks in the 10–14-year band.

4.1.3 Summary, biennial series

We conclude from the 80-year biennial records that the 11-year solar cycle can neither be confidently identified nor excluded in these data. In these records, the inherited ¹⁰Be flux represents the majority of the total ¹⁰Be flux and is much higher than the ¹⁰Be fallout flux. Expected variability in ¹⁰Be concentrations associated with 11-year-band solar variability is small in relation to measurement uncertainties, which in turn makes it difficult to recognize such variability in ¹⁰Be fallout flux estimates derived from our measurements. Spectral analysis did not display significant power in the 11-year band in fallout estimates, but if, as implied by the linear regressions discussed above, true variability in ¹⁰Be fallout is smaller than measurement errors in reconstructed ¹⁰Be fallout, we would not expect to detect it even if it was present. Therefore, our strategy of using the 11-year solar cycle to test the hypothesis that ¹⁰Be fallout variability is recorded in NAVC sediment did not yield a significant result: we have not conclusively proven that 11-year variability is present, but we also showed that, because of the poor signal-to-noise ratio induced by relatively high concentrations of inherited ¹⁰Be, we cannot prove that it is absent.

4.2 Decadal series

We now describe results from sampling two ∼ 1000-year varve sequences at decadal resolution. The purpose of these sample sets is to determine whether or not the centennial-scale variability in ¹⁰Be fallout fluxes evident in ice core records can be identified in NAVC sediments and, if identified, matched to ice core ¹⁰Be records. During some of the time periods represented by these sequences, in particular near 14 000 years ago, ice core ¹⁰Be records show centennial-period variations in ¹⁰Be deposition with a relative magnitude exceeding 50 %, substantially greater than expected from the 11-year solar cycle. Thus, it is more likely that these longer-period, larger-amplitude variations could be distinguished from measurement noise in NAVC sediments.

Figure 16 Sedimentological data and Be concentrations for decadal-resolution records from adjacent, overlapping sequences of glacial varves at Glastonbury (GL) and Kelsey–Ferguson (KF).

4.2.1 Decadal series of glacial varves from Glastonbury and Kelsey–Ferguson sites

Overlapping sequences of glacial varves from the Glastonbury (GL) and Kelsey–Ferguson (KF) sites in north-central Connecticut (Figs. 1, 2; Table 1; Fig. 16) span NAVC years 3658–4675, approximately corresponding to 17 200–16 200 yr BP. This set of samples has 20–30-year spacing; each individual sample includes 10–15 contiguous varves, and samples collected in the range of overlap of the two sections sample nearly exactly the same varves, although some pairs of samples are offset by 1–2 years. KF is 8 km north of GL and therefore more ice-proximal, but varves at KF are also slightly thinner because the site is shallower and therefore subject to a lower sedimentation rate. The interval sampled was deposited while the glacier was at least 30 km to the north of KF. The mean grain size of matching samples from the same set of varves is coarser at KF, which is consistent with the site being shallower and more ice-proximal (Fig. 16). ⁹Be concentrations are not significantly different between paired samples, but ¹⁰Be concentrations are systematically slightly higher at GL, resulting in ¹⁰Be∕⁹Be ratios that are systematically higher at GL than KF in the period of overlap between the two sample sets.

Figure 17 Be fluxes calculated for KF and GL decadal records. Upper panels: total ⁹Be and ¹⁰Be fluxes to sediment. Third panel: ¹⁰Be fallout flux estimates using Eq. (1) for constant R_S. Lower panel: combined estimate computed by linearly scaling the estimate from KF to force agreement in the period of overlap of the two records.

Figure 17 shows estimates of ¹⁰Be fallout flux for these sites. When we calculated the ¹⁰Be fallout flux for the 80-year records in the previous section, we assumed that R_S was constant and variations in the ¹⁰Be fallout flux were symmetrically distributed around some mean, which enabled us to estimate R_S from linear regression of ⁹Be and ¹⁰Be fluxes. However, for a 1000-year record, these assumptions are not likely to be true. In particular, changes in R_S and f are likely to occur on timescales of hundreds of years as the ice margin recedes from a site and the sediment grain size, the size of the contributing watershed, and the balance of glacial and paraglacial sediment all change. In fact, regressions of ⁹Be and ¹⁰Be fluxes for the long records at both KF and GL yield unphysical negative intercepts, which is consistent with the hypothesis that R_S is not constant throughout this record.

Thus, lacking any other constraints on R_S, we estimated ¹⁰Be fallout fluxes using the value of R_S that we inferred from linear regression of ⁹Be and ¹⁰Be fluxes for the biennial data at KF as discussed above in Sect. 4.1 (Figs. 17 and 13). The resulting estimates (Fig. 17) have different relative magnitudes, but similar short-term variability, with common excursions around NAVC 3850–3900 and 4000–4100. The fact that reconstructed ¹⁰Be fallout variations are similar at both sites could be evidence that we are, in fact, correctly reconstructing true fallout variations. On the other hand, reconstructed fallout fluxes at these sites are weakly correlated with total fluxes, which may indicate that our procedure has not entirely separated environmental and fallout effects on ¹⁰Be fluxes. With this value of R_S, during the period of overlap of the two records in which flux estimates show limited variability between NAVC 3850 and 4000, the mean estimated fallout flux at GL is near 1×10⁷ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, whereas at KF it is near 3.5×10⁷ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Given an approximate average fallout flux for temperate mid-latitudes of 1.5×10⁶ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Graly et al., 2011,  and references therein), these imply f≃6 for GL and f≃20 for KF, which are within the range observed for lakes elsewhere. However, as we discussed in Sect. 4.1 with respect to the biennial data, the estimate of R_S used to estimate ¹⁰Be fallout fluxes affects the mean value and the relative magnitude of reconstructed fallout variations as well as the degree of correlation between total Be fluxes and reconstructed ¹⁰Be fallout fluxes, although not the pattern of variability. Thus, although our reconstructed fallout fluxes at these sites imply physically reasonable values for f, it is difficult to conclusively prove from the data themselves that they reflect true fallout variations to the lake system. A systematic difference between fallout fluxes at the two sites could be an artifact of the calculation if R_S is not, in fact, the same at both sites, although the fact that the sediment at both sites was sourced from the same ice margin makes this unlikely. Alternatively, it could represent a real difference in the delivery rate of fallout ¹⁰Be to the sediment. For example, as GL is farther from the ice margin, the sediment deposited there could have a longer residence time in the lake, which may allow for more complete scavenging of fallout ¹⁰Be. If this hypothesis is correct, it could imply that (i) the ¹⁰Be fallout flux to the sediment Q_10,Af may be quite variable throughout the lake at any given time, and also (ii) variations in sediment transport within the lake could cause variations in the fallout flux at a given site that were unrelated to the actual atmospheric fallout rate. To summarize, our estimates of ¹⁰Be fallout fluxes at these sites appear physically reasonable and show reproducible variability in excess of measurement uncertainty, but it is difficult to conclusively prove from the measurements alone that they represent true variations in ¹⁰Be fallout to the lake system.

We explored this further by comparing ¹⁰Be fallout flux estimates at these sites to coeval ice core records. For the purposes of comparison, we assume that the difference between the two records is most likely explained by differences in f between the two sites and therefore obtain an estimate of centennial-scale variability represented in both records by applying a constant scaling factor to the GL data so as to align the two records where they overlap between NAVC 3810 and 4040, then averaging the two records within the period of overlap. This procedure loses information on the relative magnitude of fallout variability, and our combined estimate has higher internal variability (∼ 50 %; see the bottom panel in Fig. 17) than is observed in ice core ¹⁰Be flux records for this time period (∼ 15 %). However, we argued above in Sect. 4.1 that the relative magnitude of our fallout flux estimates is likely unreliable due to its dependence on R_S, so it could not be used as a point of comparison with ice core records in any case. Multitaper spectral analysis of this combined record shows a significant (at 95 % confidence) spectral peak at 91 years, which matches a characteristic period observed in solar variability records (the 88-year “Gleissberg cycle”; e.g., Feynman and Fougere, 1984). This, like the similarity of fallout flux estimates from the two records, suggests that variability in our fallout flux estimate is not solely the result of measurement noise and that true centennial-scale variability in ¹⁰Be fallout may be recorded.

Figure 18 Comparison of the combined ¹⁰Be fallout flux estimate from the KF and GL sites (Fig. 16) with ice core records of ¹⁰Be fallout. (a) Scaled and centered ice core ¹⁰Be flux records from Adolphi et al. (2018) (and references therein), offset to NAVC years with a NAVC–GICC05 yr BP offset of 20 645 years (see text). “Stack” is the global ice core stack from Adolphi et al. (2018). (b) Scaled and centered combined KF–GL fallout flux estimate compared to the ice core stack with the same 20 645-year offset. (c) Same as (b), only for an offset of 21 110 years. In panels (a)–(c), because the primary plotted units of standard deviations from the mean obscure the true relative variability of each record, the additional axes at right labeled $Q/\overline{Q}$, which have line colors and weights matching each data series, show units of relative variation with respect to the mean of each data series. This highlights the fact that the apparent magnitude of variations in our reconstructed NAVC fallout flux is greater than expected from coeval variations in ice cores, although, as discussed in the text, this difference could be an artifact of the calculation procedure and we therefore used the scaled records for correlation estimates. (d) Results of correlation analysis comparing the scaled and centered KF–GL fallout flux estimate with the ice core stack. Red highlights values of the offset for which the hypothesis that the two records are uncorrelated can be rejected at 95 % confidence. Vertical lines highlight offsets inferred from radiocarbon dating of NAVC varves (20 820 or 20 645 with the GICC05–calendar year offset of Adolphi et al., 2018, for this time period) and event correlation of varve thickness with Greenland ice core temperature (20825 and 20750;  Ridge et al., 2012). All of these yield an inverse correlation between the ¹⁰Be records.

Figure 18 compares our combined ¹⁰Be fallout flux record to ice core ¹⁰Be fallout flux records compiled by Adolphi et al. (2018). Again, as we have argued that the magnitude of variability in our ¹⁰Be fallout flux estimate is unreliable, we compare scaled and centered records to focus on the pattern of variability rather than the magnitude. The overall structure of the ice core records in this time period displays periodic centennial-band variability as well as a broader low between ∼ 16 500 and 16 800 GICC05 yr BP. If we start with the NAVC–calendar year BP offset of 20 820 years inferred from ¹⁴C data (Ridge et al., 2012; Thompson et al., 2017) and accept the estimate of Adolphi et al. (2018) that the GICC05 chronology is 175 years too young in this time interval, then we expect a NAVC–GICC05 yr BP offset of 20 $\mathrm{820}-\mathrm{175}=\mathrm{20}$ 645. Aligning the composite ¹⁰Be fallout estimate from KF and GL with the ice core ¹⁰Be stack (Fig. 18b) with a NAVC–GICC05 yr BP offset of either 20 820 or 20 645 does not yield a good match: a peak in the ¹⁰Be fallout flux near NAVC 4100 aligns with a period of relatively low flux in the ice core stack, resulting in an inverse correlation between the two records. Ridge et al. (2012), based on correlation between NAVC varve thickness records and Greenland ice core temperatures, suggested NAVC–GICC05 yr BP offsets of 20 825 for events near 16 100 yr BP and 20 750 for events near 17 400 yr BP, and neither of these yields a good match either.

We searched for a better match by looking for the closest alternative value of the offset for which the KF–GL combined fallout estimate is positively correlated with the ice core stack at high confidence (Fig. 18c, d). This occurs for NAVC–GICC05 yr BP offsets in the range 20 965–21 165 years, with the best correlation at 21 110 years. This aligns a series of fallout peaks at NAVC 3700–4100 in the KF–GL record with a broad high around 17 000–17 300 GICC05 yr BP in the ice core records and a period of low variability at NAVC 4200–4600 with a broad low around 16 500–16 800 GICC yr BP. However, the correlation is significant across a 200-year range in the offset, which shows that the correlation analysis is mainly aligning the longer-period structure of the two records rather than the centennial-scale variability. Apparent correlation of ∼ 1000-year-period variations in a record that is only 1000 years long could easily be misleading if the assumption of constant R_S and f created a spurious systematic drift.

Overall, our ¹⁰Be fallout flux estimate from the KF–GL decadal series (i) shows variations that are reproducible between two sites, (ii) shows evidence for periodic variability in a band characteristic of solar variability, and (iii) shows acceptable correlation with the global ice core ¹⁰Be fallout stack using an NAVC–GICC05 offset that is close to the expected value. These observations may suggest that our fallout flux estimate is, in fact, recording global atmospheric ¹⁰Be production and fallout variations. On the other hand, the apparent correlation between the two records mostly reflects alignment of long-period variability, which is unlikely to be robust for a relatively short record and in addition could arise from incorrectly assuming constant R_S and f. Thus, although this match between NAVC and ice core fallout records is possible and is consistent with the hypothesis that our fallout flux estimate represents true variability in ¹⁰Be fallout, it does not provide strong support for the hypothesis.

4.2.2 Decadal series of glacial and paraglacial varves from Newbury

Our decadal-resolution record from Newbury comprises a sequence of glacial varves transitioning to paraglacial varves between NAVC 6630 and 8190 (Fig. 19). This set of samples has 20–30-year spacing, and each individual sample includes 10–20 contiguous varves. Ridge et al. (2012) placed the boundary between glacial and paraglacial varves at NAVC 7100, where relatively thick varves with relatively high interannual variability (high MAR) are overlain by thinner varves (lower MAR) with reduced interannual variability (Fig. 19). However, grain size and Be concentration records (Fig. 19) do not show the same sharp transition, indicating that the glacial–paraglacial transition was more likely characterized by a gradual reduction in glacial meltwater and sediment delivery beginning near NAVC 7100 and continuing for ca. 200 years. Regardless, sediments in the lower part of this section before ca. NAVC 7000 are predominantly derived from direct glacial discharge, and sediments in the upper part after ca. NAVC 7225 are predominantly derived from runoff from the deglaciated landscape. A period of low MAR beginning at NAVC 6800 is associated with the Littleton–Bethlehem Readvance, a period of ice margin expansion into the lake recognized from a moraine complex and related stratigraphic sections ∼ 40 km to the north of this site (Ridge and Toll, 1999; Ridge et al., 2012; Thompson et al., 2017). Although the distance to the ice margin presumably decreased at this time, this effect was offset by reduced meltwater production in a cooler climate that resulted in thinner, finer-grained varves.

Figure 19 Sedimentological data and Be concentrations for the decadal-resolution record at Newbury (NEWB; see Figs. 1 and 2).

Be isotope systematics are very different in glacial and paraglacial parts of the Newbury section (Figs. 20, 21). Although, as noted above, the assumptions necessary to estimate R_S from linear regression of ⁹Be and ¹⁰Be total fluxes are probably not valid for long records, regressions for glacial and paraglacial parts of the section in Fig. 20 differ in a way that is similar to that observed for the paraglacial and glacial biennial records at NHV and KF, and they are consistent with expected changes from a glacial to paraglacial environment. Glacial varves imply $R_{\mathrm{S}}=\mathrm{1.03}\times {\mathrm{10}}^{-\mathrm{9}}$ and mean(Q_10,Af) $=\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{7}}$ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (e.g., f≃15), and ¹⁰Be fluxes display high variability around the regression line, which implies that fallout ¹⁰Be represents a significant fraction of the total ¹⁰Be flux. Paraglacial varves imply $R_{\mathrm{S}}=\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{9}}$ and mean(Q_10,Af) $=\mathrm{8.4}\times {\mathrm{10}}^{\mathrm{6}}$ $\mathrm{atoms}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (e.g., f≃5), with minimal variability. This is consistent with our hypothesis, discussed above, that R_S is expected to be higher for a deglaciated landscape that may have been accumulating ¹⁰Be in soils for thousands of years than for direct glacial sediment discharge. The decrease in f between glacial and paraglacial varves at the same site is also consistent with the prediction that f should be larger if the lake watershed includes a significant area of the adjacent ice sheet, so the glacial–paraglacial transition is associated with a significant decrease in contributing area.

Figure 20 Regression of ⁹Be and ¹⁰Be total fluxes for glacial and paraglacial varves at Newbury.

Although this analysis highlights differences in ¹⁰Be systematics between glacial and paraglacial varves, application of the estimates of R_S from the regression slopes to Eq. (1) would imply unphysical negative ¹⁰Be fallout fluxes for several samples. Thus, the estimate of ¹⁰Be fallout flux in Fig. 21 uses $R_{\mathrm{S}}=\mathrm{0.9}\times {\mathrm{10}}^{-\mathrm{9}}$. Again, the value we assume for R_S affects the mean of calculated fallout fluxes and therefore the relative magnitude of the variability but does not affect the pattern of fallout variability.

Figure 21 Be fluxes calculated for the Newbury decadal record. (a) and (b), total ⁹Be and ¹⁰Be fluxes to sediment. (c), ¹⁰Be fallout flux estimate using Eq. (1) for constant R_S.

The most striking aspect of ¹⁰Be fallout flux estimates for Newbury is that variability in the ¹⁰Be fallout flux is strongly suppressed in the paraglacial section. This is a likely consequence of buffering of ¹⁰Be fallout by watershed processes in the absence of direct glacial meltwater input. As ¹⁰Be is known to accumulate in soils in nearly all environments (Pavich et al., 1984), in a paraglacial landscape, a substantial fraction of ¹⁰Be fallout deposited in the watershed is likely adsorbed to developing soils during infiltration and therefore cannot be transported into the lake. Pollen records from this time period show that the landscape was heavily vegetated (Thompson et al., 2017), which highlights the likelihood of rapid soil development and ¹⁰Be accumulation in soils in the watershed. On the other hand, ¹⁰Be deposited on the surface of an ice sheet with negative mass balance, as presumably was the case for a large extent of the southern margin of the rapidly retreating Laurentide Ice Sheet at this time, cannot accumulate. Fallout ¹⁰Be deposited on an ablation area of the ice sheet that drains into a glacial lake must all be flushed into the lake, either by supraglacial or subglacial drainage, during the summer melt season each year. Thus, if a large area of the ice sheet drains into the lake, fallout ¹⁰Be from a large area can be effectively concentrated and transported to the lake with a minimal time lag between fallout and delivery. We conclude from the striking contrast between glacial and paraglacial ¹⁰Be systematics at this site that glacial lake sediments likely record ¹⁰Be fallout variations with a significantly higher amplitude and shorter time lag than paraglacial sediments. Fallout variations in paraglacial lake sediments appear to be strongly suppressed by watershed buffering of ¹⁰Be fluxes.

The reconstructed fallout flux in the glacial section at Newbury displays high variability with several peaks that are substantially larger than measurement uncertainty, and it does not show any correlation with total Be flux, so this variability is not likely to be explained by environmental factors. ¹⁰Be fallout records from ice cores that are coeval with the Newbury glacial section also display high variability at centennial timescales during this period, and Fig. 22 compares the ¹⁰Be fallout flux estimate from Newbury to the ice core records. Note that because the variability in the ¹⁰Be fallout flux in the paraglacial section of the record is minimal and near the level of measurement uncertainty, this portion of the record has no effect on a correlation analysis and we do not consider it. The important feature of the ice core records in this time interval is a group of prominent peaks in ¹⁰Be fallout that occurred between 13 900 and 14 500 GICC05 yr BP, which are well in excess of the characteristic variability of the surrounding time period. If we accept the estimate of Adolphi et al. (2018) that the INTCAL13 and GICC05 timescales are equivalent in this interval, the best-fitting radiocarbon calibration of the NAVC (Thompson et al., 2017) would imply a NAVC–GICC05 yr BP offset of 20 820 years. Event matching between varve thickness records and ice core climate records in the time range covered by the Newbury glacial section by Ridge et al. (2012) suggested a NAVC–GICC05 yr BP offset of 20 825 years. However, these values fail to align peaks in the varve and ice core ¹⁰Be fallout records (Fig. 22) and, in fact, yield an inverse correlation between the two records. In contrast, a NAVC–GICC05 yr BP offset near 20 950 years aligns significant peaks at NAVC 6950 and 6700 with corresponding peaks in the ice core record at 13 975 and 14 300 yr BP (Fig. 22). An offset of 20 925 yields the best alignment of the peaks at NAVC 6950 and 13 975 GICC05 yr BP (the peak at NAVC 6950 dominates a correlation metric; see Fig. 22), but a value of 20 950 better aligns both peaks (as the Newbury record has 25-year spacing, the difference between 20 920 and 20 950 is not significant). Other matches are possible because the Newbury record displays two major peaks and the ice core records display three major peaks with several adjacent minor ones in this time period. For example, an offset near 20 675 would align the peak at NAVC 6700 with the peak at 13 975 BP.

Figure 22 Comparison of the ¹⁰Be fallout flux estimate from Newbury with ice core records of ¹⁰Be fallout. Plot elements are the same as in Fig. 18. (a) Scaled and centered ice core ¹⁰Be flux records from Adolphi et al. (2018), offset to NAVC years with the value of the offset inferred from ¹⁴C data (Thompson et al., 2017) and the assumption that INTCAL13 and GICC05 timescales are equivalent in this interval. (b) Scaled and centered Newbury fallout flux estimate compared to the ice core stack for the same offset. (c) Same as (b), only for an offset of 20 925 years that aligns peaks in varve and ice core ¹⁰Be (see text). As described in the caption to Fig. 18, the additional axes to the right of panels (a)–(c) show units of relative variation with respect to the mean of each data series. (d) Correlation between varve and ice core ¹⁰Be records for different offset values. This considers only the glacial part of the section; the minimally variable paraglacial section does not affect the conclusions. Vertical lines highlight other suggested offsets inferred from radiocarbon calibration (20 820) and event correlation (20 825; see text).

Overall, the most prominent feature of the decadal-resolution ¹⁰Be record at Newbury is the strong contrast between ¹⁰Be systematics in glacial and paraglacial varves. This, as argued above, indicates that ¹⁰Be fallout variability is likely to be enhanced in glacial varves by ¹⁰Be runoff from the ice sheet but suppressed in paraglacial varves by watershed buffering. The glacial part of the section displays significant centennial-scale variability in reconstructed ¹⁰Be fallout, which can be acceptably correlated with coeval ice core records using values for the NAVC–GICC05 yr offset that are close to those inferred from other evidence.

4.2.3 Decadal series, summary results

Although the paraglacial section at Newbury does not show significant variability in ¹⁰Be fallout flux, ¹⁰Be fallout flux estimates from the glacial sections at Newbury, Kelsey–Ferguson (KF), and Glastonbury (GL) display variability in excess of measurement uncertainty that can be correlated with global ice core ¹⁰Be fallout records using NAVC–GICC05 yr BP offsets that are similar to independently estimated values. Despite several unverified assumptions in our reconstruction of ¹⁰Be fallout fluxes from Eq. (1) and significant uncertainties in relating the apparent magnitudes of ¹⁰Be fallout variability to those expected from global records, these observations are consistent with the hypothesis that apparent ¹⁰Be flux variations in NAVC sediments are, in fact, recording true fallout variations.

The NAVC–GICC05 yr BP offsets that produce the best matches at Newbury and KF–GL are different, in the range 20 900–20 950 for Newbury and 20 965–21 165 at KF–GL. The match with Newbury is substantially more likely to be valid because it is based on aligning anomalously large, short-duration peaks in both records that are expected to be insensitive to possible longer-term variations in R_S or f, whereas the apparent match for the KF and GL sections is not well defined and could be a result of assumptions in the fallout calculation. In general, estimates of the NAVC–GICC05 offset obtained for different time ranges are not required to be the same because the GICC05 timescale is not exactly linear in calendar years (Adolphi et al., 2018). If the apparent matches at Newbury and KF–GL were both valid, the difference would imply that GICC05 is overcounted by ∼ 150–200 years between 14 000 and 17 000 BP. However, other estimates of the GICC05–calendar year offset (Adolphi et al., 2018,  and references therein) show the opposite, that it is undercounted by ∼ 175 years during this interval, so a GICC05–NAVC offset of 21 000 for KF–GL should be equivalent to an offset of 21 175 at Newbury, or, conversely, an offset of 20 925 at Newbury would be equivalent to 20 770 at KF–GL. Thus, nonlinearity in the GICC05 timescale does not resolve the difference in apparent matches. Another theoretical possibility that would allow both matches to be valid would be the presence of a significant time lag between ¹⁰Be fallout variations onto the land surface and deposition in NAVC sediments, with this time lag being 200–400 years longer at KF–GL than at Newbury. However, a watershed residence time and associated fallout–deposition time lag of hundreds of years would not be consistent with the observation that 50–100-year peaks in ¹⁰Be deposition at Newbury are recorded at a similar resolution as those recorded in ice core records. Overall, these arguments tend to indicate that the apparent NAVC–ice core matches at Newbury and KF–GL cannot both be valid, and, as discussed above, other arguments indicate that of these two, the apparent match at KF–GL is more likely to be spurious and the one at Newbury is more likely to be valid. In the next sections we discuss this in light of previous independent approaches to calibrating the NAVC timescale.

4.2.4 Comparison of NAVC–GICC05 yr BP offsets inferred from ¹⁰Be flux matching with NAVC–cal yr BP calibration from radiocarbon dates

In previous work, Ridge et al. (2012) and Thompson et al. (2017) estimated the NAVC–calendar year BP offset by choosing the value of the offset that minimizes the weighted mean square of residuals between observed ¹⁴C ages of plant macrofossils in individual varves and the INTCAL13 calibration curve, yielding an estimated NAVC–cal yr BP offset of 20 820 years. This procedure assumes that all ¹⁴C ages are equally likely to be accurate, so it acts to choose a value for the offset such that the mean of the residuals is zero. However, Ridge et al. (2012) also observed that the residuals from this procedure were not normally distributed and instead were characterized by a young-skewed distribution, so even though the mean of the residuals was zero, the mode was positive (see Fig. 10 in Ridge et al., 2012, and Fig. 23). They proposed that the skewed distribution was the result of small amounts of postdepositional contamination of samples by younger carbon. This would imply that the correct value of the NAVC–cal yr BP offset should be chosen so as to align the mode of the residuals, rather than their mean, with zero, which would increase the value of the offset by 100–200 years. As shown in Fig. 23, if we accept that the GICC05 timescale is equal to true calendar years near 13 500 yr BP (Adolphi et al., 2018), then the NAVC–GICC05 yr BP offset of 20 925 years obtained from ¹⁰Be flux matching at Newbury is consistent with this reasoning: ¹⁴C residuals relative to INTCAL13 with this offset have a mode equal to zero and a young tail. In other words, with this offset, the INTCAL13 calibration curve is closely aligned with the oldest radiocarbon ages observed in clusters of data at several levels in the NAVC (Fig. 23; with the exception of four old outliers that are likely the result of contamination; see Ridge et al., 2012). Thus, the NAVC–GICC05 yr BP offset of 20 925 years inferred from ¹⁰Be flux matching at Newbury is consistent with the ¹⁴C data set and in fact may better agree with the ¹⁴C data than the estimate of the offset previously obtained by aligning the mean of the ¹⁴C residuals with zero.

Figure 23 (a) Alignment of radiocarbon ages of plant macrofossils in NAVC varves (Ridge et al., 2012; Thompson et al., 2017) with the INTCAL13 radiocarbon calibration curve for values of the NAVC–GICC05 yr BP offset estimated from ¹⁰Be flux matching (Figs. 18, 22) and the assumption that GICC05 years are equivalent to calendar years. The bands outlined in red and blue are the 68 % uncertainty bounds for the INTCAL13 calibration assuming offsets of 20 925 and 21 110, respectively. The histograms in (b) and (c) show residuals defined as the difference between the measured radiocarbon age of a sample from a particular varve year and the age predicted for that year by the INTCAL13 calibration and the selected offset.

On the other hand, applying the NAVC–GICC05 yr BP offset of 21 110 obtained from ¹⁰Be flux matching to the KF–GL records results in a poor alignment between INTCAL13 and the NAVC radiocarbon data set. Predicted ¹⁴C ages with this offset are nearly all older than observed, resulting in both the mode and the mean of the residuals being less than zero (Fig. 23) and requiring that essentially all ¹⁴C ages be spuriously young due to contamination. Adding the proposed ∼ 175-year undercount in the GICC05 timescale between ∼ 13 000 and ∼ 17 000 yr BP (Adolphi et al., 2018) worsens the agreement.

We conclude that (i) the apparent NAVC–ice core ¹⁰Be flux match at Newbury is consistent with the ¹⁴C data set, but (ii) the apparent match to the KF–GL record results in a value of the offset that is too old to be consistent with the ¹⁴C data.

4.2.5 Comparison of NAVC–GICC05 yr BP offsets inferred from ¹⁰Be flux matching with those inferred from correlation of varve and ice core climate records

In addition to calibrating the NAVC to calendar years using ¹⁴C ages, Ridge et al. (2012) and Thompson et al. (2017) proposed NAVC–GICC05 offsets based on correlation between abrupt climate changes in Greenland recorded in ice core δ¹⁸O records and ice-marginal events recorded by the NAVC. The premise of this approach is that rapid climate changes likely affected central Greenland and the margin of the LIS in the northeastern US similarly. Thus, cooling in Greenland should be associated with decreases in meltwater production, which are recorded by decreased varve thicknesses, and decreased ice margin retreat rates, or potentially stillstands or readvances, at the LIS margin. Warmings in Greenland should be associated with increases in meltwater production, varve thickness, and ice margin retreat rates. In this section, we evaluate whether these previously proposed correlations are consistent with the apparent NAVC–GICC05 matches based on comparing NAVC and ice core ¹⁰Be fallout discussed in the previous sections.

Events near 14 000 yr BP. For the period near 14 000 cal yr BP, Thompson et al. (2017) showed that the NAVC–cal yr BP offset of 20 825 years derived from fitting ¹⁴C data to INTCAL13 also yielded a close correlation between Greenland ice core climate records and events in the NAVC (Fig. 24a). This correlation associates a decrease in varve thickness near NAVC 6800 and increase near NAVC 6920 (Fig. 24) at Newbury with a rapid cooling and warming in Greenland at the boundaries of stadial GI-1d at 14 025 and 13 905 GICC05 yr BP, respectively (Lowe et al., 2008; Rasmussen et al., 2014), and places the beginning of the Littleton–Bethlehem Readvance, an ice margin readvance stratigraphically bracketed by NAVC varves 6836–6987, near the beginning of GI-1d.

Figure 24 Correlation between NAVC and GISP2 climate proxy records near 14 000 GICC05 yr BP. This duplicates Fig. 14 of Thompson et al. (2017), except that the x axis is reversed for consistency with the rest of the figures in this paper. Records at the bottom are 9-year running averages of varve thickness records from Newbury (NEWB) and Perry Hill Basin (PHS; see Ridge et al., 2012). The GISP2 δ¹⁸O record (Stuiver et al., 1995; Stuiver and Grootes, 2000) is a three-sample running average on the GICC05 chronology. The gray band shows the minimum duration of the Littleton–Bethlehem Readvance inferred from varves directly underlying and overlying the readvance till (see Thompson et al., 2017). Panel (a) uses the offset proposed as the best event correlation by Ridge et al. (2012) and Thompson et al. (2017), which is the same as the offset derived from fitting to INTCAL13 by Thompson et al. (2017). Panel (b) uses the offset that provides the best match between the ¹⁰Be flux at Newbury and the ice core stack (Fig. 22). Panel (c) uses the offset that provides the best match between the ¹⁰Be flux at KF–GL and the ice core stack (Fig. 18). Note that the anomalously thick varves near NAVC 6900 (highlighted in a) are related to a lake outburst flood and therefore have no climate significance (see Ridge et al., 2012). The Greenland climate event boundaries are from Rasmussen et al. (2014).

The NAVC–GICC05 offset of 20 925 years inferred from ¹⁰Be flux matching at Newbury (Fig. 22) implies that the Littleton–Bethlehem Readvance began in a precursory cold event prior to the start of GI-1d and ended at the GI-1d termination, which is stratigraphically permissible (Fig. 24b). On the other hand, this offset does not align appropriate varve thickness changes with the boundaries of GI-1d, so it is inconsistent with the hypothesis that cold events in Greenland should be associated with decreased meltwater production at the LIS margin. Therefore, accepting a NAVC–GICC05 yr BP offset of 20 925 years would require revising our understanding of the significance of the varve thickness changes at Newbury.

The NAVC–GICC05 yr BP offset of 21 110 inferred from ¹⁰Be flux matching at KF–GL (Fig. 18) does not yield any apparent correlation between ice core and NAVC climate proxy records in this time range, instead placing the Littleton–Bethlehem Readvance in a transitional phase of GI-1e (Fig. 24c) and yielding no obvious alternative correspondence between ice core climate and varve thickness records.

Figure 25 Correlation between NAVC and GISP2 climate proxy records near 16 000 GICC05 yr BP. This duplicates Fig. 14a of Ridge et al. (2012), except that the x axis is reversed for consistency with the rest of the figures in this paper. The plot elements are the same as in Fig. 24. The gray band shows the Hatfield event, a period of slowed ice recession or a brief stillstand associated with thinner varves between NAVC ∼ 4725 and 4825.

Events near 16 000 yr BP. Ridge et al. (2012, their Fig. 14a) proposed that a NAVC–GICC05 yr BP offset of 20 825 years also aligned ice core climate and NAVC events at the time of the Hatfield event, an ice margin stillstand between NAVC 4725 and 4825. This event is associated with an abruptly bounded period of thinner varves at South Hadley (Fig. 25), and a 20 820 offset aligns this period with an unnamed minor cold phase in the ice core record (Fig. 25a). Although the NAVC–GICC05 yr BP offset of 21 110 instead aligns this event with a warm phase (Fig. 25b), a value of 21025, which is within the fairly broad range of acceptable matches at KF–GL, would align it with a different cold phase (Fig. 25c). It is also true, however, that the offset near 20 925 years inferred from ¹⁰Be flux matching at Newbury, with the proposed ∼ 150-year undercount in GICC05 in the 13.5–16 ka range (Adolphi et al., 2018), would result in an offset near the originally proposed value of 20820. Thus, event correlation around the Hatfield event is inconclusive as to whether the ¹⁰Be flux matches at either Newbury or KF–GL are valid.

Figure 26 Correlation between NAVC and GISP2 climate proxy records near 17 500 GICC05 yr BP. This duplicates Fig. 14b of Ridge et al. (2012), except that the x axis is reversed for consistency with the rest of the figures in this paper. The plot elements are the same as in Fig. 24. Varve thickness records shown are from Glastonbury (GL), Kelsey–Ferguson (KF), and South Windsor, CT (PET; see Ridge et al., 2012). The gray band shows the period of the Chicopee Readvance, which is associated with a period of thin varves at NAVC 3320–3530.

Events near 17 000 yr BP. Ridge et al. (2012, their Fig. 14a) proposed a NAVC–GICC05 yr BP offset of 20 750 years to align ice core and NAVC climate proxy records at the time of the Chicopee Readvance, which is associated with a period of thin varves near NAVC 3320–3530. They noted that this value of the offset aligns century-scale peaks in varve thickness with corresponding ice core δ¹⁸O variations during an unnamed cold phase (Fig. 26a). The value of 20 825 proposed as the best offset for event correlation around the Littleton–Bethlehem and Hatfield events fails to align these peaks and does not correlate the readvance with the cold phase (Fig. 26b). The NAVC–GICC05 yr BP offset of 21 110 years inferred from ¹⁰Be flux matching in this time range at KF–GL fails to align the Chicopee Readvance with any cold phase (Fig. 26c). However, if we accept the NAVC–GICC05 yr BP match at 20 925 from ¹⁰Be flux matching at Newbury and apply the estimate of Adolphi et al. (2018) wherein GICC05 is undercounted by ∼ 175 years between ∼ 13 500 and ∼ 17 000 BP, then a NAVC–GICC05 yr BP offset of 20 925 at Newbury implies a NAVC–GICC05 yr BP offset near 20 750 at the time of the Chicopee Readvance, which is the same as that proposed for event correlation by Ridge et al. (2012). Thus, the ¹⁰Be flux match at Newbury is consistent with event correlation around the Chicopee Readvance, but the ¹⁰Be flux match at KF–GL is not.

To summarize, the NAVC–GICC05 offset of 20 925 years suggested by correlation between NAVC and ice core ¹⁰Be records for the glacial section at Newbury may be consistent with previously proposed correlations between NAVC and ice core climate proxy records, but the NAVC–GICC05 offset near 21 110 suggested for the KF–GL section is not.

5 Conclusions

5.1 Systematics of Be deposition in varved glacial lake sediments

We found that the majority of ¹⁰Be in glacial and paraglacial varves in the NAVC is inherited and already adsorbed to sediment when supplied to the lake. This means that variation in the total ¹⁰Be flux to the lake predominantly reflects variations in the supply of inherited ¹⁰Be, which are mostly the result of grain size variations, and not variations in fallout ¹⁰Be. By using ⁹Be as a normalizing isotope, we were able to identify variations in the ¹⁰Be flux that are not explained by variations in the inherited Be supply to the lake and which we therefore attribute to variations in atmospheric ¹⁰Be fallout. In most cases, this yielded reconstructed ¹⁰Be fallout fluxes that were weakly correlated or uncorrelated with total Be fluxes, which supports the argument that this procedure is effective at separating environmental processes from fallout variations. This is an important element of this study: if we had not measured ⁹Be concentrations, it might be difficult or impossible to separate variations in inherited and fallout ¹⁰Be fluxes. Having quantified the relationship between grain size and Be concentrations in this study, however, it might be possible in future studies of the NAVC to extract variability in ¹⁰Be fallout without ⁹Be measurements by instead identifying variations in the total ¹⁰Be flux not explained by grain size, which would be similar to the approach of multivariate correlation with various sedimentological proxies applied in other studies (e.g.,  Czymzik et al., 2016). Regardless, using ⁹Be as a normalizing isotope appears simpler.

5.2 Variations in ¹⁰Be fallout flux

Our initial strategy in this study was to test the hypothesis that ¹⁰Be fallout variations are recorded in NAVC sediments by trying to identify the characteristic 11-year solar variability period in ¹⁰Be fallout variations. This strategy was not successful because of the high ratio of inherited ¹⁰Be to fallout ¹⁰Be in NAVC sediments. Expected fallout variability of 15 %–20 % around the mean associated with the 11-year solar cycle translates into variability in the total ¹⁰Be flux that has a similar magnitude as the measurement uncertainty. Thus, we were not able to distinguish short-period solar forcing of ¹⁰Be fallout variations, if present, from measurement noise.

On the other hand, in 1000-year varve sequences sampled with biennial resolution at two sites, we did observe centennial-scale variability in the reconstructed ¹⁰Be fallout flux resembling ¹⁰Be fallout variations observed in ice cores. This leads us to believe that ¹⁰Be concentrations in NAVC sediments do, in fact, record at least some prominent high-amplitude variations in ¹⁰Be fallout that can potentially be used to correlate NAVC and ice core records.

5.3 ¹⁰Be systematics in glacial and paraglacial varve sequences

An unexpected result of this study is that glacial varves are likely better suited to reconstructing ¹⁰Be fallout than paraglacial varves. Initially, we supposed that paraglacial varves would be more likely to yield a record of ¹⁰Be fallout simply because they are typically an order of magnitude thinner than glacial varves, so ¹⁰Be fallout would be less diluted by higher sedimentation rates and associated larger inherited Be fluxes. This supposition was incorrect, and we found that variability in ¹⁰Be fallout was evident in glacial varve sequences but nearly entirely suppressed in paraglacial varves.

This result is most likely explained by the different roles of the watersheds that contribute to glacial and paraglacial lakes. In a paraglacial environment in which the lake is not receiving direct glacial meltwater, the lake watershed predominantly consists of a deglaciated landscape mantled by fresh glacial sediment. A significant fraction of ¹⁰Be fallout to the watershed is therefore adsorbed to sediment in developing soils, which limits the delivery of fallout ¹⁰Be to the lake and buffers variations in fallout rates.

In contrast, the contributing watershed of a proglacial lake in direct contact with the ice margin most likely includes large areas of the ablation zone of the ice sheet. Fallout ¹⁰Be cannot be stored on the ice in the ablation zone, so each year's fallout must be flushed into the lake by supraglacial or englacial discharge during the summer melt season. We hypothesize that this phenomenon can result in substantial focusing of fallout ¹⁰Be in glacial lake sediment and a minimal lag time between fallout variations and recording of those variations in glaciolacustrine sediment records. If true, this implies that varved glacial lake sediments, if deposited in a lake with a large supraglacial watershed, may be very effective recorders of atmospheric ¹⁰Be fallout. This hypothesis could easily be tested in modern proglacial lakes with large contributing areas on the ice surface, perhaps in west Greenland.

5.4 Correlation of the NAVC to ice core records by ¹⁰Be flux matching

We found values of the NAVC–GICC05 yr offset that resulted in acceptable correlations between reconstructed variations in ¹⁰Be fallout flux in NAVC sediments and an ice core ¹⁰Be flux stack for two varve sequences: a sequence at Newbury deposited near 14 000 yr BP and overlapping sequences at Kelsey–Ferguson (KF) and Glastonbury (GL) deposited ca. 16 000–17 000 yr BP. The NAVC–GICC05 yr BP offsets derived from these matches are not the same for both sections and also differ by 50–200 years from previously proposed values of the NAVC–calendar year offset inferred from radiocarbon dating of NAVC sediments and event matching between NAVC varve thickness records and Greenland ice core climate records.

Matching ice core ¹⁰Be records with ¹⁰Be fallout flux variations inferred from the older KF and GL varve sequences implies a NAVC–GICC05 yr BP offset in the range 20 965–21 165 years, with the best correlation at 21 110 years. However, the match between the two records mainly relies on aligning millennial-scale trends in the two records rather than distinct peaks associated with centennial-scale variations. As the KF–GL record is only 1000 years long, this could lead to spurious results. In addition, centennial-scale variations in the ¹⁰Be fallout flux in ice core records during this time period are relatively subtle. This match is not consistent with the NAVC–cal yr BP offset inferred from fitting NAVC ¹⁴C data to the INTCAL13 calibration curve, even when likely differences between the GICC05 ice core timescale and the true calendar year timescale are taken into account. Although offsets in the range implied by this match yield an acceptable correlation between ice core climate records and NAVC varve thickness records for events near NAVC 4800 and 16 000 GICC05 yr BP (Fig. 25), they are not consistent with the expected relationship between NAVC and ice core climate proxy records for younger or older periods (Figs. 24, 26). Overall, the apparent match between the KF–GL ¹⁰Be fallout flux and ice core ¹⁰Be does not appear consistent with most independent evidence and is unlikely to be valid.

Matching the ¹⁰Be fallout record from the glacial section at Newbury with ice core ¹⁰Be records, in contrast, relies on correlation of prominent large-magnitude, centennial-scale variations in both records and implies a NAVC–GICC05 offset in the range 20 900–20 940 (Fig. 22), with the best correlation at 20 925 years. This value of the offset is consistent with the ¹⁴C data set and, if the likely undercounting of ice core years in the GICC05 timescale between ∼ 13 000 and 17 000 yr BP is considered, also consistent with correlations between NAVC and Greenland ice core climate proxy records for periods surrounding the Hatfield event near NAVC 4800 and 16 000 GICC05 yr BP (Fig. 25) and the Chicopee Readvance near NAVC 3400 and 17 400 GICC05 yr BP (Fig. 26). This value of the offset yields an acceptable correlation between the stratigraphic boundaries of the Littleton–Bethlehem Readvance and the GI-1d stadial (Fig. 24) but does not clearly align the associated period of decreased varve thicknesses, which is interpreted to reflect decreased meltwater production during a cold period, with stadial boundaries. Thus, the 20 925-year NAVC–GICC05 yr BP offset inferred from ¹⁰Be fallout matching at Newbury is consistent with the majority of independent evidence and may be valid, but a reassessment of the climate significance of varve thickness variations associated with the Littleton–Bethlehem Readvance would be required to resolve all conflicts with independent data.

5.5 Weaknesses and strategy for improvement

The weakest part of our argument that the correlation between the ¹⁰Be fallout flux in NAVC sediments at Newbury and the ice core stack is valid is that the correlation is based on a short record comprising only 600 glacial varves that display two prominent peaks in reconstructed ¹⁰Be fallout flux. The ice core stack, on the other hand, includes three major fallout peaks between 14 000 and 15 000 GICC05 yr BP and a number of relatively high-amplitude secondary features. In addition, we only observed these ¹⁰Be fallout peaks at one section, so we cannot prove by replication between sections that they represent true fallout variations. Finally, the 25-year sampling resolution at Newbury is long relative to the size of the expected fallout peaks. For example, one of the most prominent apparent peaks is only represented by a single sample, and this effect alone could potentially explain much of the apparent ∼ 100-year difference between our proposed correlation based on the ¹⁰Be flux and that based on correlation of climate records in previous work. Another significant weakness derives from the property of our linear model for reconstructing the ¹⁰Be fallout flux that the apparent amplitude of reconstructed fallout variations displays a strong dependence on the value chosen for the inherited ¹⁰Be∕⁹Be ratio R_S, so the amplitude of reconstructed fallout variations cannot be used as a point of comparison with ice core records. These weaknesses could in part be addressed by (i) higher-resolution sampling at the Newbury section to better resolve apparent fallout peaks and (ii) replication of the ¹⁰Be fallout record from additional varve sequences that overlap with Newbury in the range NAVC 6000–7000. Several overlapping sections exist, including the Perry Hill Basin site shown in Fig. 24, that could be used for this purpose. Higher-resolution ¹⁰Be measurements from older parts of the NAVC might be helpful in resolving whether or not the NAVC–ice core ¹⁰Be match at KF–GL is or is not valid, but they are less likely to be successful because ¹⁰Be fallout variability in the ice core stack between 15 000 and 18 000 GICC05 yr BP has a lower amplitude than that in the 13 000–15 000 GICC05 yr BP range. Overall, collecting higher-resolution ¹⁰Be data on replicate sections of glacial varves from the portion of the NAVC between approximately NAVC 5800 and the glacial–paraglacial transition near NAVC 7000 is likely the best strategy for resolving the remaining ambiguous results of this study.