Spatially averaged stratigraphic data to inform watershed sediment routing: An example from the Mid-Atlantic United States

Individual storms can mobilize, transport, and deposit large volumes of fluvial sediment over a few days, but these sediment-transporting events may only occur every few years or decades. Sediment deposited during such storms may remain in alluvial deposits for decades, centuries, and even millennia before being remobilized and carried farther downstream by subsequent storm events. Thus, sediment transport through watersheds encompasses multiple time scales, from short-term events to episodic transport and storage over thousands of years.

While geologists recognize the range of time scales that control sediment transport and delivery (Brierly and Fryirs, 2005; Martin and Church, 2004; Meade, 2007), most sedimentrouting models currently used for watershed management focus on transport by individual storm events and rely on observations of sediment flux from stream gauging stations for calibration (Borah et al., 2008; Shenk and Linker, 2013). These models generally do not explicitly account for sediment storage processes, ignoring the loss of sediment through deposition, the decadal to millennial time scales associated with alluvial storage, and the gain of sediment to stream channels through erosion and remobilization of stored sediment. The use of these models, and their reliance on the short-term record of stream gauging for calibration, reflects a conceptualization of sediment delivery as a short-term, event-driven process that is inconsistent with geological data, which demonstrate the importance of centennial to millennial time scales. Models that neglect geological time scales will necessarily fail to capture transient processes associated with climate and land-use changes, watershed restoration strategies, and other important drivers.

Advances in both modeling and data analysis will be needed to effectively include longer-term sediment storage processes in predictions of watershed-scale sediment transport. Fortunately, substantial progress in model development has already occurred. A few process-based models that account for sediment storage have been developed (Carroll et al., 2004; Coulthard and Macklin, 2003; Delft 3D-FLOW, 2020; Nicholas, 2013), but these are difficult to apply over spatial scales of 10³ km and temporal scales of >10³ yr because of model complexity and limitations of computing power and data storage.

Generalized models that represent alluvial deposits as spatially averaged “storage reservoirs” are more promising. Storage reservoirs are defined regions in space that receive sediment through deposition and release sediment through erosion (Dietrich et al., 1982). They can be considered “spatially averaged” because storage reservoirs typically occupy regions larger than the scale over which specific erosional and depositional processes operate. Sediment contained within storage reservoirs can be characterized by mass, grain-size distribution, age distribution (Bradley and Tucker, 2013), carbon content (Torres et al., 2017, 2020), concentration of heavy metals (Pizzuto, 2020), and other characteristics. The time during which sediments remain in storage before being remobilized can be estimated as well; this is termed the transit or storage time distribution (Bolin and Rodhe, 1973). While storage reservoirs play important roles in hydrology, ecology, and other disciplines, geology is perhaps uniquely positioned to exploit concepts of reservoir theory because alluvial storage reservoirs are readily observed and sampled, and the disciplines of stratigraphy and geochronology provide well-established methods with which to define the characteristics of alluvial sediment reservoirs and their temporal variability.

Spatially averaged storage reservoirs can be utilized in sediment-routing models in different ways. One approach is to represent an entire river corridor by a series of connected reservoirs, with sediment transfers between them governed by rules extracted from sediment budgets, with all stored sediment equally subject to erosion (Kelsey et al., 1987; Malmon et al., 2002, 2003). Another approach is to couple a traditional process-based model for channel sediment routing to off-channel storage reservoirs. Deposition into sediment reservoirs is related to sediment concentration and other variables extracted from the channel transport model, while erosion from sediment reservoirs is determined by sediment age and other characteristics (Lauer, 2012; Lauer and Parker, 2008a, 2008b; Pizzuto, 2020).

A current limitation of reservoir-based sediment routing was identified by Lauer and Parker (2008a, p. 10), who used an exponential distribution to represent age and storage time distributions in their model, noting that “we know of no other distributions having been applied to river floodplains.” A few subsequent field studies concluded that storage time distributions for river corridors are likely to be heavy-tailed distributions, with a greater probability of remobilizing younger sediment than older sediment (Konrad, 2012; Lancaster and Casebeer, 2007; Lancaster et al., 2010; Miller and Friedman, 2009 [as cited by Bradley and Tucker, 2013]; Torres et al., 2017, 2020). However, field data available to characterize the age distribution of alluvial deposits over geologic time scales are limited.

Full development of sediment-routing models that account for sediment storage will require site-specific field data for model conceptualization and calibration. Data are needed to quantify the characteristics of sediment storage reservoirs, which must include estimates of the stored sediment mass and its age and storage time distributions, all of which can vary spatially and through time. While these data are needed by modelers, only geologists trained in stratigraphy and geochronology have the ability to create these data sets. At present, however, stratigraphic studies of geologic history are typically not designed to generate quantitative results useful for testing and developing mass-conserving numerical models.

This paper presents a case study that illustrates how geologic data can be organized to facilitate development and testing of reservoirbased sediment-routing models. New field data and previously published studies of river corridor sedimentation associated with European settlement of the Mid-Atlantic region were combined to document temporal changes in stored sediment mass and age distributions from the Holocene to the present. These results provide a useful perspective on the geologic consequences of the well-documented watershed disturbances in the region associated with anthropogenic activity, but, more importantly, they also provide a methodological template for characterizing time-dependent alluvial storage to aid development and calibration of long-term sedimentrouting models. Thus, while sediment routing is intended to be the ultimate application of the results of this study, the primary focus here is on presenting the geologic data needed to quantify sediment storage.

Relatively few studies have documented the details of Holocene, pre–European settlement valley-bottom sedimentation in the Mid-Atlantic region, though continental-scale assessments suggest that sedimentation rates were low before European colonization (Kemp et al., 2020; Wilkinson and McElroy, 2007). Jacobson and Coleman (1986) interpreted presettlement stream deposits as representing sinuous channels with low-lying floodplains composed of thin, fine overbank deposits underlain by thin, laterally accreted sands and gravels. Walter and Merritts (2008, p. 299) described valley bottoms composed of “small anabranching channels within extensive vegetated wetlands that accumulated little sediment.” Wegmann et al. (2012) offered a similar interpretation, citing presettlement gravels overlain by thin (<0.5 m) deposits of organic-rich clay-to-silt loam that they interpreted as wetland or waterlogged soil layers. Several authors described areas where presettlement valley bottom sedimentation had been enhanced by activities of indigenous peoples (Merritts et al., 2011), though James (2019) inferred that these areas were limited spatially. Presettlement valley bottoms were probably heavily impacted by beavers (Brush, 2009; Butler and Malanson, 2005; Ruedemann and Schoonmaker, 1938).

Overlying the presettlement deposits, there is a thick layer of sand, silt, and clay. These deposits, often termed “legacy sediments,” are typically associated with European colonization, the expansion of agriculture, and widespread mill damming in the Mid-Atlantic region (Costa, 1975; Jacobson and Coleman, 1986; Pizzuto, 1987; Walter and Merritts, 2008) and elsewhere in the United States (Happ et al., 1940; Kemp et al., 2020; Knox, 1987, 2006; Trimble, 1971, 1983; Wilkinson and McElroy, 2007). Legacy sediments are variably interpreted as resulting from overbank deposition arising from the increased flooding and erosion associated with poor agricultural practices (Jacobson and Coleman, 1986), and sedimentation in millponds (Dow et al., 2020; Merritts et al., 2011, 2013; Wegmann et al., 2012; Walter and Merritts, 2008). Contemporary erosion of legacy sediments is frequently cited as a potential waterquality issue for Mid-Atlantic watersheds (Merritts et al., 2013; Miller et al., 2019; Noe et al., 2020).

While initial studies of Mid-Atlantic valley bottom sedimentation did not emphasize contemporary vertical accretion on floodplains, recently published measurements indicate that overbank processes continue at significant rates. Hupp et al. (2013), Schenk et al. (2012), and Noe et al. (2020) measured ongoing sedimentation at many sites in the Mid-Atlantic region using clay pads, while Bain and Brush (2005) and Pizzuto et al. (2016) documented contemporary vertical accretion on valley bottoms using fallout radionuclides. Others inferred contemporary sediment storage on floodplains from sediment budgets (Donovan et al., 2015; Smith and Wilcock, 2015). These observations of vertically accreted deposits are supplemented by many studies that document ongoing deposition of coarse-grained, laterally accreted floodplains in the region (Donovan et al., 2015; Jacobson and Coleman, 1986; Pizzuto et al., 2014; Walter and Merritts, 2008).

Despite the extensive research summarized above, important questions regarding U.S. Mid-Atlantic valley bottom deposits remain unanswered. The relative proportions of pre- and postsettlement deposits, while documented in many localized case studies, remain poorly constrained on a regional basis, which may have important implications for water quality. Ongoing contemporary vertical accretion, while widely documented, has not been presented in a stratigraphic context, so the relative importance of contemporary floodplain sedimentation is not widely recognized. Documenting sediment storage time scales associated with these deposits is important to understand the routing of sediment through watersheds (Pizzuto, 2014), but few studies have attempted to quantify the duration of storage and to incorporate the resulting temporal lags into models, even though such models play a central role in watershed management in the region (STAC, 2005).

Following Lauer and Parker (2008a, 2008b), Lauer (2012), and others, Pizzuto (2020) described a modeling framework that includes channelized suspended sediment transport and an off-channel floodplain storage reservoir. Sediment routing downstream is accomplished by traditional mechanistic equations for suspended sediment transport, with deposition into storage related to sediment supply (e.g., suspended sediment concentration), sediment deposition velocity, the scale of the floodplain reservoir, and additional variables that account for changes in the accommodation space available for deposition through time. A separate mass balance equation accounts for stored sediment mass through time and its age distribution. The following key equation specifies removal of stored sediment through time:

where dER(T,t) is the erosion rate (mass/time) of stored sediment in an age category T at time t, dM(T,t) is the portion of the total stored mass of sediment of a given age T at time t, r(t) is an erosion rate constant (1/time) that may vary with time, and e_p(T) is a dimensionless “area exposure function” (hypothesized by Pizzuto [2020] to be time-invariant to reflect the architecture of sediment storage) that defines the exposure of stored sediment of age T to active erosion. By specifying erosion based on exposure of stored sediment by age, Equation 1 sets the time scale associated with sediment storage, and it also determines the storage time distribution of stored sediment. Because erosion is related to the stored mass dM by age category, the age distribution must be known as a function of time to use Equation 1. This information is provided by a mass balance equation formulated for the storage reservoir. The model thus predicts the stored sediment mass and its age and storage time distributions, all of which vary with time.

Figure 1 illustrates aspects of this model that might be suitable for routing sediments in a typical Mid-Atlantic watershed before and after European settlement. The watershed itself is represented by a spatial domain of river corridors. Because sediment storage time scales are “long,” computations should begin in the Holocene, before European settlement of the region; in Figure 1, a “time 0” of 1000 yr in the past is somewhat arbitrarily selected. Water discharge and sediment supplied to the upstream ends of the network from hillslope erosion and other processes represent the required boundary conditions; these must be specified over the entire time period of model predictions. To initiate model computations of erosion using Equation 1, the stored sediment mass and its age distribution must be specified as initial conditions. Then, the model steps forward through time, computing channelized transport, rates of deposition and erosion, stored sediment mass and its age distribution, and sediment yield at the downstream end of the watershed.

The modeling framework of Figure 1 requires geologic data to initiate model computations, and these data can also play an essential role in model calibration. Data are needed to define the mass of stored sediment and its age distribution 1000 yr ago, both of which are needed to initiate erosion rate computations. Once model predictions are available, data can also be used to calibrate parameters that control transport, erosion, and deposition. Typically, direct observations of erosion and deposition are unavailable for past conditions, and it is only the stored sediment that can provide useful information for calibration. For present conditions, a sediment budget can be constructed (Reid and Dunne, 2016) that can be used to evaluate the model’s ability to predict current rates of transport, erosion, deposition, and sediment yield. Once the model has been calibrated using data covering the past 1000 yr (including contemporary observations), in principle, it should be capable of evaluating sediment movement into the future (Fig. 1), though boundary conditions resulting from future land-use and climate changes may be difficult to specify.

The focus of this study was on assembling the necessary data to implement this modeling framework, which simulates millennial-scale sediment routing. Of course, sediment transport and storage processes are not only controlled by long-term processes of sediment storage but are also influenced by shorter-term processes related to individual hydrologic events, anthropogenic forcing, and other processes. These shorter-term controls are addressed by many studies (e.g., Borah et al., 2008) and will not be discussed further here.

This study combined new results with previously published data. New data were collected at second-to fifth-order stream reaches of southeastern Pennsylvania, Maryland, and Virginia (Table 1; Fig. 2). Most study reaches were located within the Piedmont Physiographic Province, a region of moderate relief underlain by Paleozoic metamorphic rocks with a humid temperate climate (Noe et al., 2020). One study site along the South River, Virginia, was located within the Valley and Ridge Physiographic Province (Bingham, 1991). Land uses included a mix of agriculture, suburban development, forest, and pasture. All of the study watersheds have been affected by a past history of watershed disturbance associated with European colonization and settlement likely beginning in the seventeenth and eighteenth centuries (James, 2019; Noe et al., 2020). Rivers in this region are gravel-bedded streams with banks composed of sand, silt, and clay. Fluvial processes are typically impacted by localized valley confinement and exposed bedrock (Pizzuto et al., 2018; Bodek et al., 2021).

Field sites were located at the White Clay Creek in southeastern Pennsylvania, and along Difficult Run and the South River in Virginia (Fig. 2). At the White Clay Creek site, we sampled a vertical sequence of deposits exposed on an eroding bank (site 1) to fully document the depositional history of these sediments, and we also determined contemporary floodplain accretion rates at five additional sites (Table 1). Detailed descriptions of these sites can be found in Bodek (2020) and McCarthy (2018). We also collected data at two sites along Difficult Run in northern Virginia, one near U.S. Geological Survey (USGS) stream gauging station Difficult Run above Fox Lake near Fairfax, Virginia (USGS #01645704), where the drainage basin area is 14 km², and the other just downstream of Leesburg Pike, where the drainage basin area is 117 km². Gellis et al. (2017), Hupp et al. (2013), and Schenk et al. (2012) provided detailed descriptions of these sites in Difficult Run. Finally, we collected data at two locations along the South River, denoted RRkm 4.75 and RRkm 5.58, where the terminology indicates the distance measured along the channel of the South River in 2005 downstream (north) from a footbridge across the South River to an industrial plant in Waynesboro, Virginia. These two sites have been featured in previous studies by O’Neal and Pizzuto (2010), Pizzuto (2012), Pizzuto et al. (2010, 2016, 2018), and Washburn et al. (2018), while additional useful background information regarding the South River has been provided by Pizzuto and O’Neal (2009) and Rhoades et al. (2009).

In addition to new data obtained for this study, stratigraphic and geochronologic data were also extracted from published data sources (Table 2) from study areas in Delaware, Maryland, Pennsylvania, and Virginia (Fig. 2). These data included interpretations of the thicknesses of pre- and postsettlement deposits, ages of presettlement deposits, and depths to chronostratigraphic horizons representing sediments deposited in 1950 and 1880. Site locations and study methodologies are described in the original publications.

At our field sites, we described stratigraphic sections and obtained samples for sediment dating and supplemented these studies with geomorphic mapping to document sediments resulting from recent lateral accretion and colonial mill damming. We combined published sediment dating with our new data to create a spatially averaged vertical accumulation history of typical valley-fill deposits for the Mid-Atlantic region. These data, when appropriately transformed, yielded estimates of sediment age distributions at key times in the past. Finally, to emphasize the implications of these data for sediment routing and delivery, illustrative sediment travel-time distributions were computed using a previously published mathematical model, a naïve assumption of a balanced steady-state sediment budget, and typical sediment budget parameters for our study area.

To estimate rates of sediment accumulation, we sampled sediments from eroding streambanks at the South River and Difficult Run and from eroding streambanks and floodplain cores at the White Clay Creek site. Sediments were described in the field according to depth below the floodplain surface, grain size, organic matter content, color, presence or absence of stratification, and other characteristics. Samples were obtained for ¹⁴C dating, optically stimulated luminescence (OSL) dating, and analyses of the fallout radionuclides (FRN) ²¹⁰Pb and ¹³⁷Cs.

We also created a geomorphic map at the Miller Heights location in Difficult Run, Virginia (Fig. 2). We classified the valley bottom into two surficial landform categories following Allmendinger et al. (2007), Jacobson and Coleman (1986), and others. The first category represents a spatially extensive, higher-elevation valley bottom surface comprising most of the valley flat, while the second category is a slightly lower surface typically found on the insides of actively migrating meander bends. According to Allmendinger et al. (2007), Jacobson and Coleman (1986), and others, the more extensive, higher-elevation surface formed by vertical accretion from the Holocene to the present, while the lower-elevation surface formed by lateral accretion processes during the twentieth and twenty-first centuries (Donovan et al., 2015). Hupp et al. (2013) indicated that a colonial-age mill dam was once located just upstream of the USGS stream gauge at the Miller Heights site. We did not observe any evidence of this structure in the field, but we located the dam on our geomorphic map based on results presented by Hupp et al. (2013).

Sediments were dated using a variety of methods. Presettlement deposits were dated using ¹⁴C and OSL dating techniques. Younger deposits were dated using ¹³⁷Cs, ²¹⁰Pb, and dendrochronology. Sedimentation rates were extracted from the ¹³⁷Cs and ²¹⁰Pb data using a numerical model described by Pizzuto et al. (2016). All dating methods were consistent with standard procedures, which are discussed in detail in the Supplemental Material.¹

We used Undatable software (Lougheed and Obrochta, 2019a, 2019b) to produce age-depth curves for each of our study sites. Undatable uses a Bayesian framework to account for uncertainty in age dating and sampling depth to produce age-depth curves with well-quantified uncertainties. Radiocarbon dates were entered into Undatable in calendar years, so Undatable’s ¹⁴C calibration routines were not needed. Depths (with uncertainties) associated with the year 1950 A.D. were obtained from modeling of ²¹⁰Pb profiles. Depths associated with peak ¹³⁷Cs activities were obtained from a numerical model of sediment and radionuclide accumulation (Pizzuto et al., 2016), corrected for vertical displacement, and were assigned (with uncertainties) to the year 1963 A.D., the year of peak atmospheric ¹³⁷Cs activity associated with nuclear bomb testing (He and Walling, 1996). Age-depth bootstrapping was enabled for all points except for those at the top and bottom of each stratigraphic section (Lougheed and Obrochta, 2019b).

From the age-depth curves, we extracted the thicknesses and associated uncertainties of the following chronostratigraphic horizons: 1950–present, 1750–1950, and all sediments older than 1750. We term these chronostratigraphic units “modern” (1950–present), “legacy” (1750–1950), and “presettlement” (older than 1750) deposits. These chronostratigraphic horizons are somewhat similar to the presettlement, agricultural, and very recent stratigraphic units defined by Jacobson and Coleman (1986); they are intended to identify sediments deposited before, during, and after watershed impacts associated with European colonization of the Mid-Atlantic region.

The specific dates defining the boundaries of these units are somewhat arbitrary, because the timing of the watershed changes they are meant to capture is difficult to determine precisely. The year 1750 is a reasonable approximate date for the initiation of colonial-age watershed disturbance, but settlement patterns, dam construction, and deforestation and agricultural development likely were not coeval throughout the entire Mid-Atlantic region. Jacobson and Coleman (1986) and Walter and Merritts (2008) noted that streams began to form floodplains by lateral accretion sometime early in the twentieth century. The timing of this transition is poorly known. The year 1950 is selected to represent the boundary between legacy and modern sedimentation largely based on analytical convenience: Within the framework of Undatable, this chronostratigraphic horizon can be readily determined using ¹³⁷Cs and ²¹⁰Pb data.

We combined the results obtained from our own field studies with published stratigraphic data for the Mid-Atlantic region. Our goal was to better document the thicknesses of presettlement, legacy, and modern chronostratigraphic units, and to more clearly define the maximum ages of presettlement deposits. However, very few previous studies provided sufficient dated sections to define all three chronostratigraphic units; rather, they more typically defined one or two of the units, often without any age control (e.g., descriptions of the thicknesses of pre- and postsettlement deposits at a location).

To circumvent the problem of incomplete data, we developed a method to synthesize complete stratigraphic sequences from the available data (Table 2). The first step in the approach was to summarize stratigraphic data as empirical frequency distributions. We combined our own new data with published observations to create frequency distributions of presettlement thickness, postsettlement thickness, depth to the 1950 chronostratigraphic horizon, and depth to an 1880 A.D. chronostratigraphic horizon defined by the history of chromium mining in Maryland from Bain and Brush (2005). The contact between presettlement and postsettlement deposits was assumed to represent a mean calendar date of 1750 A.D., with a standard deviation of 25 yr to account for regional variations in the timing of watershed disturbance associated with colonial settlement. Available dates from presettlement deposits were used to define the distribution of possible maximum ages for this stratigraphic unit.

Synthetic stratigraphic sections were created by randomly sampling empirical cumulative frequency distributions and then assembling stratigraphic columns from the data using an algorithm designed to ensure that stratigraphic superposition was always satisfied (Fig. 3). Empirical frequency distributions were randomly sampled by (1) selecting a uniform random variate between 0 and 1, and (2) obtaining the value of the variable (age or depth) by interpolation from the empirical frequency distribution (Press et al., 1992). Stratigraphic sections are represented by the age-depth coordinates of the following four horizons: the base of the section (often, but not always, presettlement in age), the contact between presettlement and postsettlement deposits, and the 1880 and 1950 chronostratigraphic horizons. Because few data were available to define the depth of the 1880 chronostratigraphic horizon, it was only included in 20% of the synthetic stratigraphic sections. Sedimentation rates for presettlement, legacy, and modern deposits were also computed from the synthetic agedepth data.

Stratigraphic sections created by this approach can be viewed as a statistical summary of regional stratigraphic data, assembled without regard to geography, watershed position, geologic setting, or any other spatial attribute. The potential limitations that arise through this extreme spatial averaging are addressed in the Discussion.

According to the conceptual framework of Figure 1, the age distribution of stored sediment is a vital component of sediment routing where storage is significant. The cumulative age distribution of stored sediment can be obtained from the synthetic age-depth sections by simply converting the depth coordinate to a fractional depth by dividing a given depth value by the maximum depth of each section. This simple conversion creates a cumulative frequency distribution of age for each synthetic stratigraphic section.

The synthetic stratigraphic sections represent the history of vertically accreted sediment in Mid-Atlantic river valleys, typically beginning before European colonization and continuing to the present. These vertically accreted deposits represent a variety of depositional environments and processes, including wetlands (Jacobson and Coleman, 1986; Walter and Merritts, 2008), overbank deposition on floodplains (Jacobson and Coleman, 1986), as well as vertical accretion imposed by mill dams in areas outside of the locations of present river channels (Walter and Merritts, 2008; Wegmann et al., 2012). Identifying specific depositional environments is beyond the scope of this study, but it is particularly important to recognize that these stratigraphic sequences represent a continuous (though episodic) record of vertical accretion without significant erosion.

Because these stratigraphic sections represent a record of continuous vertical accretion, it is straightforward to estimate age distributions of stored sediment for any time in the past. This can be achieved by selecting a specific time at which to determine the past age distribution, removing all the stored sediment deposited after that time (e.g., all sediment younger than the specified age), and rescaling the depth axis from 0 to 1, with 0 representing the sediment surface at the specified age and 1 representing the base of the section. To illustrate how age distributions of vertically accreted sediment changed following European colonization, we present age distributions for 1000 yr ago and 1900 A.D. in addition to contemporary age distributions that represent all preserved deposits. Depths associated with 1000 yr ago and 1900 A.D. were obtained from the synthetic stratigraphic sequences by linear interpolation.

The synthetic vertical sections (and new sections created from our field data) do not entirely account for some important deposits that have been identified in Mid-Atlantic river corridors. These include filling of preexisting channels upstream of mill dams constructed soon after European colonization, the subsequent removal of these deposits early in the twentieth century when most of the mill dams failed (Merritts et al., 2011; Pizzuto and O’Neal, 2009; Walter and Merritts, 2008), and floodplains formed by lateral accretion beginning around 1900 (Allmendinger et al., 2005; Jacobson and Coleman, 1986).

A specific example from our field site near Miller Heights in Difficult Run (Figs. 2 and 4) illustrates how age distributions can be corrected for these additional processes. This site is useful because a geomorphic map and other data are available to specify the mass of sediment stored upstream of a mill dam and sediment created by lateral migration since 1900. The resulting corrected age distributions are valid only for this location; a general approach for including these deposits must await more widespread geomorphic mapping. However, the results are useful to demonstrate how such corrections can be made.

To adjust the age distribution for 1900 A.D. to account for sediment deposited in the channel upstream of a mill dam, the age of the dam and accumulation of sediment upstream through time must be specified. This information is not known in Difficult Run in detail; rather, some reasonable assumptions are made here to illustrate both the approach and the magnitude of the resulting corrections to the age distribution. The dam is assumed to have been constructed in 1800 A.D. Accumulation of sediment upstream is assumed to have filled the channel through time exponentially, such that the fraction of the bankfull channel filled with sediments followed the equation 1 – e^–t/25, where t is the number of years since 1800. These assumptions provide the data required to estimate the mass of stored sediment in the channel in any age category, and these masses of relatively young sediment can then be added to the masses of sediment stored by vertical accretion to obtain the corrected age distribution for sediment stored in 1900.

A procedure was also designed to correct the present age distribution for laterally accreted sediment since 1900, with the additional assumption that all the mill dam deposits stored in the channel had been completely removed at the time of our sampling in 2017 (as indicated by our geomorphic mapping). However, preliminary computations indicated that the volume of laterally accreted sediments was insignificant compared to the volume of vertically accreted deposits in our study area (Fig. 4), and therefore corrections for lateral accreted deposits are unnecessary. Similar conclusions were reached by Pizzuto et al. (2014) based on data from the Good Hope Tributary in Maryland. Donovan et al. (2015) reported average lateral migration rates of ~0.025 channel widths per year in Maryland, further suggesting that lateral reworking has been minimal in the twentieth and twenty-first centuries.

While the age distributions created using methods described above are useful, storage time distributions are more closely related to sediment delivery, as they define how long particles remain in storage before being remobilized and carried farther downstream. Storage retards sediment delivery, creating significant lag times between watershed management actions and improvements in water quality (Meals et al., 2010; Science and Technical Advisory Committee, 2005). While rigorous evaluation of sediment travel times is beyond the scope of this paper, it is still useful to provide an approximate evaluation of the sediment travel times that are implied by the age distributions created during this study.

To assess sediment travel times, it is first necessary to estimate the storage time distribution, which typically is not identical to the age distribution of stored sediment. As noted by Lancaster et al. (2010) and Lancaster and Casebeer (2007), the storage time distribution during a particular time interval may be defined by the age distribution of sediment eroded during that time interval. This provides a useful conceptual framework for relating age and storage time distributions and also for estimating storage time distributions by dating sediments actively exposed to erosion (for example, in eroding cutbanks).

Many studies in the Mid-Atlantic region describe actively eroding near-vertical streambanks with exposures of pre- and postsettlement deposits (Donovan et al., 2015; Hupp et al., 2013; Jacobson and Coleman, 1986; Merritts et al., 2013; Schenk et al., 2012). Our geomorphic mapping in Difficult Run (presented in Results below) suggested that most eroding banks at this location are composed of vertically accreted presettlement, legacy, and modern deposits. The age distribution of these deposits was defined by the methods of the previous section, and because they are actively eroding, the contemporary storage time distribution is equivalent to their age distribution. Because these observations have been widely noted throughout the region, the equivalence between the contemporary age distribution of vertically accreted deposits and the storage time distribution is likely a general conclusion.

To illustrate sediment travel times consistent with the age distributions obtained from our analysis, we followed an idealized approach wherein watershed sediment budgets are treated as balanced, with erosion equal to deposition, and storage distributions and sediment budgets are assumed to be temporally invariant. These computations are based on Pizzuto et al.’s (2017) stochastic approach, where storage time scales are defined by a known probability distribution. In this framework, sediment moves downstream at a constant rate (consistent with the concept of a constant effective water discharge), the frequency of which is assessed using an intermittency factor (Paola et al., 1992). The downstream drift during transport is set at 0.5 m/s, and it is assumed that transport occurs 1% of the time, for a total yearly downstream transport of 158 km/yr (i.e., 0.5 m/s converted to km/yr ´ 0.01). Pizzuto et al. (2014, 2017) noted that sediment travel times are primarily controlled by storage for long transport distances, so results are insensitive to the choice of the drift velocity. For each kilometer of transport, the probability of storage is evaluated as the fraction of the sediment load stored per kilometer, q, which is specified from a contemporary sediment budget. Computations presented here assumed a value of 0.01 km^–1 for q, a typical value for the Mid-Atlantic region (Pizzuto et al., 2014; of course, available data indicate that q can vary considerably, the implications of which are explored further in the Discussion). If sediment is stored at a particular location, then the duration of storage is determined from a known probability distribution of storage times. Because it is convenient analytically, and also because it is consistent with our data, contemporary storage times were assumed to be exponentially distributed, with a mean storage time (residence time) of 300 yr. To illustrate how travel-time distributions could vary with watershed size, travel-time distributions are presented for distances from 10 to 500 km.

The mathematical basis for travel-time computations was discussed by Pizzuto et al. (2014, 2017) and is omitted here for brevity. A derivation of the relevant equations is summarized in the Supplemental Material.

Vertically accreted deposits dating from Holocene to present underlie most of the valley bottom at the Miller Heights site (Fig. 4) in the Difficult Run watershed. Difficult Run has eroded some of these deposits during the twentieth and twenty-first centuries and created new deposits through lateral accretion, but contemporary laterally accreted deposits are limited in extent (Pizzuto et al. [2014] presented a similar geomorphic map of the Good Hope Tributary in Maryland). The extent of channel-fill deposits behind the colonial-age mill dam was inferred, as these deposits had been entirely eroded away by 2017, when the mapping was completed. Eroding banks often, but not exclusively, occur along the outer banks of bends, where erosion is slowly removing vertically accreted deposits adjacent to the channel (Fig. 4).

Eroding banks sampled for dating consisted of 1.7–4.2 m of sand, silt, and clay overlying basal sand and gravel (Fig. 5). The upper meter was typically composed of massive sandy mud with a few thin layers of muddy sand or sand, but surficial deposits at the White Clay Creek site consisted primarily of sand. Colors were yellowish brown (10 YR 5/6), dark yellowish brown (10 YR 4/4), or dark reddish brown (2.5 YR3/3). Little organic matter was preserved. At South River RRkm 4.75 and the White Clay Creek site, presettlement deposits were capped by a decimeter-thick, dark-gray (10 YR 3/6, 10 YR 3/1) mud unit generally interpreted to represent a buried A horizon developed on the presettlement floodplain surface (Jacobson and Coleman, 1986). At the Leesburg Pike site in Difficult Run, a log road was preserved at a depth of ~3 m (also described by Hupp et al., 2013), indicating that these deposits postdated European colonization. At all the sites, lower deposits above basal sand and gravel consisted of mud and muddy sand with occasional layers of sand, were poorly stratified, and displayed frequent rust-colored mottling. Organic content was generally low, though twigs, leaves, and other organic materials were occasionally found (Difficult Run–Miller Heights, White Clay Creek), often immediately above basal sand and gravel layers.

New OSL and radiocarbon dates ranged in age from 17,150 yr to 130 yr (Table 3; Fig. 5). The oldest dates were obtained from banks of the South River; basal dates from the Difficult Run sites ranged from 630 yr to 190 yr, while the basal ¹⁴C date from the White Clay Creek site 1 was 616 yr. Multiple dates from individual sections nearly always obeyed superposition, except for the two lowest dates from the South River RRkm 4.75 and the two lowest dates from Difficult Run Leesburg Pike: At these two locations, the lower of the two dates was slightly younger than its immediate upper neighbor. Dates obtained from Difficult Run at Leesburg Pike ranged in age from 130 to 210 yr, consistent with the required postsettlement age of these deposits defined by their location above the buried log road. Dates below buried A horizons at the South River RRkm 4.75 and the White Clay Creek sites, all greater than 600 yr old, are consistent with the interpretation of this paleosol as representing the presettlement floodplain surface (Jacobson and Coleman, 1986). The OSL date overlying the buried soil at the South River RRkm 4.75 site, with an age of 310 yr, might initially appear too old to be postsettlement in age, as this would require the watershed disturbance associated with European settlement to predate 1700 A.D. in this area of Virginia, but Wegmann et al. (2012) cited a similar age for some basal postsettlement deposits in North Carolina. Furthermore, basal postsettlement sediments are known to be time-transgressive (James, 2019), and they have been rarely dated directly, so the 300-yr-old age obtained here is not unreasonable.

Fallout radionuclide (FRN) data from upland reference sites at Difficult Run and White Clay Creek showed activity concentrations that decreased rapidly with depth (Fig. 6). At Difficult Run, peaks in ¹³⁷Cs and ²¹⁰ Pb activity concentrations occurred just below the surface, with activity concentrations decreasing rapidly thereafter. The steady-state atmospheric deposition rate of ²¹⁰Pb at the Difficult Run upland site is 0.015 (±0.001) Bq/cm²/yr (0.90 dpm/cm²/yr, where dpm is disintegrations per minute), slightly higher than measured annual atmospheric fluxes of 0.72 dpm/cm² reported from Lewes, Delaware (Hartman, 1987), 0.78 dpm/cm² from the Chesapeake Bay (Kim et al., 2000), and 0.79 dpm/cm² (1983) and 0.85 dpm/cm² (1984) from Norfolk, Virginia (Todd et al., 1989). Data from the White Clay Creek upland site are similar to those from Difficult Run, but they display greater variability. The steady-state atmospheric deposition rate of ²¹⁰Pb estimated from the White Clay Creek data is 0.021 ± 0.0004 Bq/cm²/yr (1.3 dpm/cm²).

FRN data from the three exposed banks displayed consistent patterns with depth. Peak ¹³⁷Cs activity concentrations occurred well below the surface at depths ranging from ~20 cm at Difficult Run–Miller Heights (Fig. 7) to ~40 cm at the other sites (Figs. 7 and 8). Activity concentrations of ²¹⁰Pb_ex decreased approximately exponentially with increasing depth, with welldefined profiles at the two Difficult Run sites but extensive scatter at the White Clay Creek site. Activity concentrations of ²¹⁰Pb_ex decayed to near 0 at depths of ~40 cm at all three sites (Figs. 7 and 8).

Modeled profiles of FRN distribution with depth at all the sites defined varying model parameters and sedimentation rates. Modeling at reference sites defined rates of atmospheric deposition, bioturbation for ¹³⁷Cs and ²¹⁰Pb_ex, and downward advective movement of ¹³⁷Cs. Parameter ranges were generally similar at exposed banks, where sedimentation was included as an additional process, though at some sites, parameter ranges for the floodplain sites were expanded beyond those used to model FRN accumulation at the upland reference sites. Bioturbation coefficients ranged from 2 ´ 10^–8 to 2 ´ 10^–6 kg²/cm⁴/yr (Table 4), while advective movement of ¹³⁷Cs ranged from 0.04 to 0.16 cm/yr. Advective rates of ¹³⁷Cs would displace the ¹³⁷Cs peak 2–9 cm from 1963 to 2017, indicating the importance of accounting for this process. Mean sedimentation rates obtained from model calibration ranged from 0.29 to 1.10 cm/yr.

Sedimentation rates from dendrochronology supplemented results obtained from ¹³⁷Cs and ²¹⁰Pb_ex (Table 5; Fig. 9) for the White Clay Creek site. Sediment ranging in thickness from 8 to 44 cm overlaid basal roots of trees ranging in age from 6 to 38 yr, resulting in sedimentation rates that varied from 0.5 to 5.2 cm/yr. Neglecting the extreme outlier of 5.2 cm/yr, the mean sedimentation rate from dendrochronology is 1.0 cm/yr (standard deviation = 0.5 cm/yr), while the uncertainty resulting from propagating the individual standard deviations for each tree is 1.3 cm/yr. Obtaining representative sedimentation rates by dendrochronology is often limited by the spatial distribution of suitable trees for analysis, which could be a concern at site 1 at the White Clay Creek (Fig. 9), where sampled trees were either located along the river bank or in two small areas upstream. Nonetheless, the mean sedimentation rate from dendrochronology of 1.0 cm/yr is within the range of 0.44–1.10 cm/yr determined from the ¹³⁷Cs and ²¹⁰Pb_ex data (Table 4).

Age-depth relationships from the Undatable software illustrated the history of sedimentation through time and documented the thicknesses of presettlement, legacy (1750–1950), and “modern” (post-1950) deposits at the five study sites (Fig. 10). Four of the five age-depth curves showed an abrupt increase in slope around the time of European settlement, representing increasing sedimentation rates associated with watershed disturbance. This increase in sedimentation rate was not displayed by results from Difficult Run–Leesburg Pike, because no deposits of presettlement age were exposed at this site. At both South River sites and the White Clay Creek site, presettlement deposits exceeded legacy and modern deposits in thickness. Modern sediments deposited since 1950 were significant at all five sites and ranged in thickness from 22 ± 2 cm at Difficult Run–Miller Heights to 43 ± 10 cm at Difficult Run–Leesburg Pike.

Sedimentation rates extracted from the agedepth relationships generally showed increasing rates through time (Fig. 11). Presettlement sedimentation rates were very low at four sites where these deposits are present, while legacy and modern sedimentation rates were much higher. Unexpectedly, modern sedimentation rates were approximately equal to or exceeded legacy sedimentation rates at all the sites except Difficult Run–Leesburg Pike, where legacy sedimentation rates exceeded modern sedimentation rates by a factor of ~3. Sedimentation rates obtained from FRN analyses agreed reasonably well with those obtained using dendrochronology at Difficult Run (Schenk et al., 2012) and at the White Clay Creek site (Figs. 9 and 11; Tables 4 and 5).

Rates of vertical accretion determined from ²¹⁰Pb and ¹³⁷Cs analyses at five additional sites in the White Clay Creek watershed ranged from 0 to 0.75 cm/yr (Fig. 12; Table 1). Uncertainties at two locations (sites 8 and 9) included the possibility of 0 sedimentation, and at a third site (AJP-2), a sedimentation rate of 0 is the preferred interpretation. These sedimentation rates were combined with published values to increase the available data documenting contemporary floodplain accretion.

Published data document thicknesses and ages of presettlement deposits, the thicknesses of postsettlement deposits, and the depths to layers representing deposition in 1880 and 1950 (Fig. 13). Postsettlement deposits are the thickest of these, with a median thickness of 89 cm, considerably exceeding the median thickness of presettlement deposits of 55 cm. Approximately 7% of stratigraphic sections have no presettlement deposits preserved. The median age of presettlement deposits is 900 yr, but ages range from less than 300 yr to over 17,000 yr. The median depth of the 1950 surface is 18 cm, while the corresponding value for 1880 is 26 cm. Thickness values for all the stratigraphic units overlap one another; for example, some of the sediments deposited since 1950 exceed the thickness of some postsettlement deposits. This illustrates that the thicknesses of these stratigraphic units can vary over a wide range from place to place.

The variability of stratigraphic relationships implied by the data in Figure 13 is highlighted by the distributions of three stratigraphic contacts extracted from 5000 synthetic stratigraphic sequences displayed in age-depth coordinates in Figure 14. The base of presettlement deposits can be encountered at depths of less than a decimeter to depths of nearly 400 cm. Similar variability is displayed by the age-depth coordinates of the other two contacts. Average age-depth coordinates of these contacts, however, define a smooth upward trend that documents increasing sediment storage through time.

To summarize the relative contribution of presettlement, legacy, and modern deposits to the total thicknesses of the 5000 synthetic stratigraphic sequences of Figure 14, the thicknesses of the three stratigraphic units were expressed as fractions of the total thickness of each synthetic section. The distributions of fractional thicknesses (Fig. 15) indicate that the relative proportions of the three units vary widely. The median fraction of presettlement deposits is 0.39, but fractional thicknesses range from 0 to ~0.7. Legacy sediments typically comprise 47% (median) of the total thickness of the synthetic stratigraphic sequences, but thicknesses can range from 0 to 100%. The median fraction of modern sediments is 0.11, but the fractional thicknesses of individual sections also vary from 0 to 100%.

Distributions of sedimentation rates of the three stratigraphic units obtained from the 5000 synthetic sections are presented in Figure 16. The median presettlement sedimentation rate is 0.06 cm/yr, consistent with the results from the five field sites shown in Figure 11, but sedimentation rates can vary from 0 to >1 cm/yr. The median sedimentation rate of legacy sediments is 0.34 cm/yr, which is slightly higher than the median sedimentation rate of modern sediments, which is 0.25 cm/yr, but sedimentation rates for both units can range from 0 to >1 cm/yr.

Contemporary probability density functions of sediment age of 5000 synthetic sections are highly variable, but they display a prominent peak for stored sediment ~200 yr old, representing enhanced sedimentation rates associated with European settlement (Fig. 17A). Sediments associated with legacy sedimentation display a frequency approximately twice that of younger deposits, a trend that is most clearly shown when the probability density function is computed by pooling all 5000 sections together (solid line in Fig. 17A). An exponential age probability density function with a mean of 300 yr is consistent with some of the observed age probability density functions (dashed line in Fig. 17A), but it lacks the peak associated with legacy sediments, highlighting the importance of time-varying deposition rates in controlling the present age distribution of stored sediments in the Mid-Atlantic region.

The age probability density function reconstructed for 1900 A.D. is similar to the present age probability density function, except that the enhanced frequency associated with legacy sediments is now associated with younger sediments (Fig. 17B). This is reasonable, because in 1900 A.D., most of the legacy deposits were <100 yr old, and these sediments have aged in place another 100 yr since 1900 A.D. when analyzed for present conditions. Correcting the probability density function of “all” 5000 sections preserved in 1900 A.D. for millpond sediments nearly doubles the frequency associated with legacy sediments, indicating the potential importance of these deposits for past age probability density functions reconstructed at times when millpond deposits were in place.

The age distribution of presettlement deposits is presented as a complementary cumulative distribution function in Figure 17C, along with a complementary exponential distribution with a mean of 2160 yr included for reference. The exponential distribution represents the observed data well for the youngest ~90% of preserved deposits, but older sediments are more abundant than is predicted by an exponential distribution. This “heavy tail” results from a few relatively rare, older deposits, particularly those from basal deposits of the South River, and it is broadly consistent with heavy-tailed alluvial storage functions proposed by Bradley and Tucker (2013), Lancaster et al. (2010), Pizzuto et al. (2017), and Torres et al. (2017, 2020).

Cumulative travel-time distributions resulting from an exponential age distribution with a mean of 300 yr are presented in Figure 18. For a travel distance of 10 km, more than 90% of the sediments are delivered within 30 yr, but the percentage delivered in 30 yr drops to less than 10% at distances of ~300 km (Fig. 18B). At distances of a few hundred kilometers, the median travel time iŝ10² yr, with some sediments requiring 10³ yr to travel to the outlet of a watershed (Fig. 18A). At a distance of 500 km, the median travel time is >10³ yr.

While previous studies have often emphasized the prevalence of legacy sediments in Mid-Atlantic valley-fill deposits (Dearman and James, 2019; Walter and Merritts, 2008; Wegmann et al., 2012), our new field data and synthesis of previously published data demonstrate that presettlement and modern stratigraphic units are significant as well. The median fraction of presettlement deposits in reconstructed synthetic sections is 0.39, somewhat less than the corresponding value of 0.47 for legacy sediments. Sediments deposited since 1950 are also significant, representing a median fraction of 0.11 of the preserved deposits. These results suggest that legacy sediments (as defined here) are typically not the predominant stratigraphic unit in Mid-Atlantic floodplains, which challenges the current narrative (e.g., Walter and Merritts, 2008; and others) and has significant implications for management of these landforms (Miller et al., 2019).

Sedimentation rates for the three stratigraphic units defined in our analysis are also informative. The low rate of accumulation of presettlement deposits (median 0.06 cm/yr) is consistent with results of previous studies, as is the rapid increase in sedimentation associated with European colonization (median 0.34 cm/yr; Kemp et al., 2020; Knox, 1987, 2006; Wilkinson and McElroy, 2007). The relatively high rate of contemporary sedimentation (median 0.25 cm/yr) has also been documented previously (Hupp et al., 2013; Noe et al., 2020; Pizzuto et al., 2016; Schenk et al., 2012), but it has rarely been included in a comparative analysis of rates associated with older deposits of the Mid-Atlantic region, and these modern deposits have also not often been considered in a stratigraphic context in previous studies.

While the average thicknesses and sedimentation rates of the three stratigraphic units are important, their extensive variability is also significant. Any of the three stratigraphic units can be absent locally, and thicknesses can vary substantially. Ages of presettlement deposits, when present, can vary from a few hundred years to more than 17,000 yr. While some variability in our results could arise from neglecting coherent geographic variations in sediment storage processes (Dearman and James, 2019; Knox, 1987; Magilligan, 1985) or from imprecise dating or stratigraphic interpretations, much of the variability is likely an essential feature of river corridor sedimentation that should be explicitly accounted for in sediment routing and watershed management planning.

It is also important to consider the limitations of the available data when evaluating our results. Our study was a quantitative evaluation that attempted to represent valley corridor deposits over a large region, while the available data are not extensive (as is illustrated by Tables 1–3). For example, there are almost no stratigraphic data available for the largest rivers in the region, and dated stratigraphic sections are relatively rare compared to the thousands of kilometers of river lengths in our study area. Furthermore, the limited data that are available are not randomly distributed spatially but rather reflect the interests and biases of individual scientists. More extensive data would likely reveal additional and possibly different results from those presented here.

Contemporary and past age distributions extracted from synthetic stratigraphic sections exhibited a number of important features. First, contemporary age distributions were highly variable, suggesting that the ages of stored sediment will change considerably from section to section. Second, contemporary age distributions displayed a prominent peak representing the substantial pulse of legacy sediment deposited since 1750. While all alluvial deposits by definition reflect past conditions to some degree, the peaked contemporary age probability density function clearly illustrates how events of the last few centuries continue to dominate the sedimentary record of Mid-Atlantic river corridors.

Characteristic time scales associated with pre- and postsettlement deposits are considerably different. Representative exponential functions shown in Figures 18A and 18C, while not precisely equivalent to observed relationships, capture some important features of these differing time scales. The mean of 2160 yr for the exponential function fit to presettlement deposits suggests a long residence time for these deposits, while the mean of 300 yr, mostly fit to legacy and recent deposits, indicates a much shorter residence time. These differences in characteristic storage time scales have important consequences for watershed management now and into the future.

Travel-time distributions of Figure 18 were included to illustrate several important principles. The first is simply that contemporary age distributions of stored sediment that span centennial time scales, interpreted as storage time distributions, imply lengthy travel times even though suspended sediment is carried rapidly downstream in the channel during storm events. Many of the travel times illustrated in Figure 18 greatly exceed time scales envisioned by watershed restoration plans, here represented by the 30 yr travel time of Figure 18B. Results of Figure 18B imply that watershed restoration activities will only benefit areas within a few tens of kilometers downstream within 30 yr, reinforcing the need to include the methods of this study in sediment-routing models designed to evaluate restoration scenarios. A second implication of these results, emphasized in previous studies (e.g., Pizzuto et al., 2014, 2017), is that over long transport distances (e.g., hundreds of kilometers), the timing of sediment delivery is essentially independent of transport processes within the channel (e.g., individual storms), and it entirely depends on the time sediments spend in storage between transport events. Finally, it is also noteworthy that the results of Figure 18 were obtained analytically; a complex numerical model was not required.

Despite the useful illustrative principles derived from these travel-time distributions, it should be remembered that they are presented only as idealized examples. They are only strictly valid for a single sediment budget scenario that may or may not broadly represent the Mid-Atlantic region, and they were created using a very simplistic exponential approximation to the many storage time distributions extracted from our synthetic stratigraphic sequences. Furthermore, only a single value for the fraction of sediment exchanged per kilometer, q, was used, even though data from the region indicate that q can vary by at least an order of magnitude (Pizzuto et al., 2014). Pizzuto et al. (2014) demonstrated that mean travel times are proportional to q (their equations 1 and 22), so travel-time predictions are clearly sensitive to local variations in sediment budgets, a factor that is not explicitly reflected in the results of Figure 18. Finally, the computations assume a steady-state system where erosion and deposition are balanced and unchanging through time. None of these conditions can be a reasonable approximation of contemporary sediment transport processes in the Mid-Atlantic region.

The results presented herein highlight some of the key geologic data required to account for sediment storage in a watershed-scale sediment-routing model. The erosion equation presented as Equation 1 requires the mass of stored sediment and its age distribution to be specified at the beginning of model computations, which occurs 1000 yr ago for the framework of Figure 1. The methods adopted here defined both quantities through the results of Figures 14 and 17. As model computations proceed through time, the accumulation of stored sediment and age distributions are predicted, allowing data for subsequent time periods (e.g., age distributions for 1900 and the present) to be used for model evaluation and calibration.

All of our results exhibited substantial variability. As noted above, most variability likely arises from natural fluctuations in sedimentary processes and is therefore irreducible. Because variability is expected and unavoidable, geoscientists and modelers should develop formal processes to assess how uncertainty in initial conditions and data for calibration affects uncertainty in model predictions.

While our results illustrate some essential geologic data required by models, other information will also be needed to implement the framework of Figure 1. Sediment supplied to the channel network by hillslope erosion and water discharge must be determined through time so that boundary conditions are fully specified. A variety of anthropogenic influences should be accounted for, including the dates of mill dam construction and demise, their locations within each river corridor, and their influence on sedimentation through time. Spatial variations in relevant parameters and processes, neglected here, are likely important factors in specific watersheds and must also be assessed.

Spatially averaged stratigraphic and geochronologic data, summarized as age-depth relationships, provided a quantitative summary of vertically accreted valley bottom sedimentation from Holocene to present in the Mid-Atlantic region. Stratigraphic data were divided into three chronostratigraphic units: sediments deposited before widespread watershed disturbances associated with European colonization, termed presettlement deposits, sediments deposited as a result of watershed disturbances between 1750 and 1950, termed legacy sediments, and modern sediments deposited after 1950. While the thicknesses of all three stratigraphic units are characterized by exceptional variability, the median fractional thicknesses of presettlement, legacy, and modern deposits are 0.39, 0.47, and 0.11, indicating that legacy sediments are not necessarily the predominant stratigraphic unit in Mid-Atlantic valley-fill deposits. Presettlement deposits accumulated at a median rate of 0.06 cm/yr, while sedimentation rates accelerated dramatically to 0.34 cm/yr following European settlement. Modern sedimentation rates are only slightly lower than those associated with legacy sediments, with a median value of 0.25 cm/yr. The sedimentation rates of all units are highly variable: Trends in sedimentation rates defined by median values may not be apparent in many individual sections.

Reconstructed age-depth relationships were recast as age distributions, and, when combined with geomorphic mapping and assessment of patterns of erosion through time, as storage time distributions as well. Age distributions for the presettlement period, in 1900 A.D., and at present are highly variable, but they illustrate systematic variations through time. The presettlement distribution reflects preserved deposits as old as 18,000 yr and a slow reworking time scale (residence time) of ~2160 yr. By 1900, rapid deposition of legacy sediment created an age distribution dominated by sediments younger than 150 yr in age (as indicated by the peak in the age distribution for these deposits depicted in Fig. 17B), though deposits from the early Holocene continue to be represented. The present age distribution also exhibits increased frequency of legacy sediments, reflecting the continued preservation of these deposits in Mid-Atlantic river corridors. The residence time for contemporary sediments is ~300 yr, indicating a transition from a presettlement fluvial system that slowly reworked stored sediment to a more vigorous contemporary fluvial system with a considerably shorter residence time.

Travel-time computations were presented to illustrate potential relationships between the age distribution of stored sediments and the distribution of sediment travel times through large Mid-Atlantic watersheds. Inferred travel times ranged from a few years to more than 1000 yr, suggesting that storage can create long lag times between watershed restoration to reduce sediment loading and delivery of benefits to watershed outlets.

An explicit representation of sediment storage is clearly needed to improve watershed-scale sediment-routing models. Such models will require geologic data for development and calibration. The methods and results described here illustrate how stratigraphic and geochronologic data can be analyzed and organized to achieve these goals.

Financial support was provided by the Delaware Watershed Research Fund award DWRF-16-109, National Science Foundation grant EAR-1424969, and Petroleum Research Fund award PRF 57428-ND8. K. Skalak and A. Benthem assisted with stratigraphic description, geochronology, and analysis of samples from the South River and Difficult Run. S. Mahan performed OSL dating of samples from the South River and Difficult Run. M. Sherif and A. Pearson contributed 137Cs and 210Pb analyses of samples from White Clay Creek. J.E. Pizzuto participated in all aspects of the study. Stephanie Stotts assisted with dendrochronology at site 1 of White Clay Creek, Neil Sturchio assisted with fallout radionuclide analyses, and Joey George created the geomorphic map of Difficult Run. Perceptive and helpful reviews of a previous version of this manuscript were provided by Michael Church and an anonymous reviewer; the authors are grateful for their contributions. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government. OSL data availability: OSL data sets related to this article have not yet been approved by the USGS. When approved, the data will be archived in USGS Science Base and can be accessed by typing in the name of the associated article. Interested researchers are encouraged to access the data at https://www.usgs.gov/products/data/data-releases.