Floodplain Sediment Storage Timescales of the Laterally Confined Meandering Powder River, USA

As sediment is transported through river corridors, it typically spends more time in storage than transport, and as a result, sediment delivery timescales are controlled by the duration of storage. Present understanding of storage timescales is largely derived from models or from field studies covering relatively short (≤102 year) time spans. Here we quantify the storage time distribution for a 17 km length of Powder River in Montana, USA by determining the age distribution of eroded sediment. Our approach integrates surveyed cross‐sections, analysis of historical aerial imagery, aerial LiDAR, geomorphic mapping, and age control provided by optically stimulated luminescence (OSL) and dendrochronology. Sediment eroded by Powder River from 1998 to 2013 ranges from a few years to ∼5,000 years in age; ages are exponentially distributed (r2 = 0.78; Anderson‐Darling p value 0.003). Eroded sediment is derived from Powder River's meander belt (∼900 m wide), which is only 1.25 times its meander wavelength, a value reflecting valley confinement rather than free meandering. The mean storage time, 824 years (95% C.I. 610–1030 years), is similar to the time required to rework deposits of Powder River's meander belt based on an average meander migration rate of ∼1 m/yr, implying that storage time distributions of confined meandering rivers can be quantified from remotely sensed estimates of meander belt width and channel migration rates. Heavy‐tailed storage time distributions, frequently cited from physical and numerical modeling studies, may be restricted to unconfined meandering rivers.

Storage times have been estimated using a variety of methods. Sediment budgets typically provide estimates of the mass of stored sediment and the erosional and depositional fluxes of sediment into and out of storage; these data provide estimates of sediment residence times, but they cannot define the full distribution of sediment storage times (Lauer & Parker, 2008b;Repasch et al., 2020). Residence times can also be inferred from travel times derived from the spatial distribution of radionuclide concentrations, but these data also do not define the full distribution of sediment storage times (Dosseto et al., 2006;Lauer & Willenbring, 2010;Repasch et al., 2020;Whiting et al., 2005). Time series of aerial imagery and dendrochronology have provided estimates of storage time distributions for floodplain sediments based on areal patterns of age and erosion. These studies are typically limited to short timescales of 10 2 years or less, and often suggest that sediment storage times are exponentially distributed with an equal probability of erosion for all ages (Everitt, 1968;Nakamura, 1986;Nakamura & Kikuchi, 1996;Nakamura et al., 1987), with some studies also proposing power-law distributions (Konrad, 2012;Miller & Friedman, 2009;Pizzuto et al., 2017).
Physical and numerical modeling studies can potentially evaluate storage over longer timescales; these methods often indicate that storage time distributions are heavy-tailed, meaning that the distribution has a non-finite first and/or second moment such that younger sediment has a higher probability of being eroded than older sediment. Heavy-tailed distributions are common in natural phenomena (Caers et al., 1999a(Caers et al., , 1999b, and frequently appear in studies of fluvial sediment dynamics. For example, Ganti et al. (2011) observed that in a flume study by Sheets et al. (2007), a heavy-tailed distribution of storage times was produced and best fit by the heavy-tailed Pareto distribution. Pareto distributions also appear to well-represent the age distribution of sediments reworked by some numerical models of river meandering (Torres et al., 2017), though Bradley and Tucker (2013) argue that the heavy-tailed Lévy distribution better represents results obtained in their simulation of river meandering based on the CHILD landscape evolution model because of the diffusive-like nature of channel position over time. Finally, Bradley (2017) directly observed heavy-tailed storage time distributions through a decade-long field tracer experiment. Results from these studies, though compelling, remain untested by field studies over longer timescales of 10 2 -10 4 years.
Determining the mathematical form of storage time distributions is an important step toward interpreting and understanding how storage influences sediment routing. If storage time distributions are exponential, for example, analysis of storage processes is greatly simplified. Exponential distributions are scaled by a single parameter that can often be readily estimated, while the heavy-tailed distributions have multiple parameters that are difficult to quantify and interpret (Bradley & Tucker, 2013). In steady-state systems characterized by an exponential storage time distribution, storage time and age distributions are identical (Bolin & Rodhe, 1973), so defining one distribution immediately determines the other. Exponential storage time distributions can also lead to robust analytical methods for sediment routing (Malmon et al., 2002), so the implications of sediment storage on sediment delivery can be assessed without complex numerical modeling.

10.1029/2021JF006313
3 of 21 Two approaches are available to define complete storage time distributions from field observations. The first approach is exemplified by Moody (2017), who derived a decadal time series of storage time distributions from comprehensive field surveys. However, due to the labor involved, such comprehensive data are rare. The second approach was used by Lancaster and Casebeer (2007) and Lancaster et al. (2010), who dated mass-wasting and other deposits exposed in eroding montane riverbanks, thereby providing a direct assessment of the age distribution of sediments just prior to erosion. Lancaster and Casebeer (2007) and Lancaster et al. (2010) proposed power laws to represent their observations. Similarly, Torres et al. (2020) dated particulate organic carbon in fluvial deposits and riverine suspended sediments using 14 C and similarly found that younger material is preferentially remobilized with a Pareto distribution providing the best fit to their data.
These studies reveal important gaps in our understanding of sediment storage times. Field data available to estimate sediment storage times over geologic timescales are limited. While physical and numerical models suggest that sediment storage times are heavy-tailed, with younger deposits being preferentially eroded in relation to older deposits, few field observations are available to verify this hypothesis. Furthermore, the mathematical form of storage time distributions should also vary systematically with geomorphic setting, with the lateral extent of valley sedimentation exerting an important control. For example, Repasch et al. (2020) document relatively short storage timescales within the meander belt of the Rio Bermejo in the Amazon foreland of Argentina, and much longer timescales associated with alluvial fan deposition outside of the river's meander belt. These observations suggest that valley confinement (Fryirs et al., 2016;Lewin & Brindle, 1977) may limit the maximum timescale of stored sediments. In confined valleys, the age limit of stored sediment could conceivably be scaled by a simple ratio of meander belt width to channel migration rate, while in unconfined fan settings, the age limit of stored sediment may be more closed tied to avulsion frequency (Repasch et al., 2020).
The purpose of this study is to fill this knowledge gap by documenting the annual to millennial floodplain storage time distribution of a largely undisturbed meandering river system. We achieve this by mapping and dating eroding floodplain deposits of the meandering Powder River in southeastern Montana. Our objectives are (a) to test the hypothesis that the sediment storage time distribution of meandering rivers is "heavy-tailed," and (b) to assess the role of valley confinement on the age distribution of stored sediment.

Study Area
Powder River is a free-flowing meandering river on the high plains of the western United States. It drains 34,700 km 2 of northeastern Wyoming and eastern Montana, flowing 600 km north-northeast before joining the Yellowstone River ( Figure 1). This study focuses on an approximate 90 km reach between the U.S. Geological Survey (USGS) Moorhead (USGS #06324500) and Broadus (USGS #06324710) stream gaging stations (Figure 1). Within this study area, we performed detailed geomorphic mapping and other observations centered within a shorter, 17 km length of Powder River's valley between Moorhead and Broadus ( Figure 1).
Powder River's hydroclimate and geomorphic character are typical of the region. The climate is semi-arid. Powder River is subject to different types of floods: ice-jam floods in early spring, snowmelt flooding in May and June, flooding due to convective summer thunderstorms, and occasionally, frontal flooding in September and October (Moody, 2019). At the Moorhead gaging stations, the mean-annual flow was 12.5 m 3 /s between 1930 and 2010 and the mean annual bankfull discharge was 170 m 3 /s (Moody, 2019) near the gaging station at Moorhead. Powder River carries a sediment load consisting of gravel, sand, silt, and clay, with an annual suspended load of 2-3 million tons at Moorhead . The low-flow channel width is approximately 50 m, the slope averages 0.001, and the sinuosity ranges from 1.1-3.6 (Martinson, 1984). Powder River has only a few manmade structures (bridges and local pump sites only in use during the irrigation season) allowing for essentially unrestricted meandering and channel migration (Gay et al., 1998). Schook et al. (2017) report a 184year record of meander migration based on cross-section surveys since 1975, aerial imagery since 1939, and tree ring analyses since 1830. Meander migration rates averaged 1.62 m/yr from 1830 to 1939, decreased to 1.52 m/yr from 1940 to 1975, and decreased further to 0.81 m/yr from 1975 to 2014, suggesting that a meander migration rate of ∼1 m/yr represents a reasonable average migration rate for Powder River during the thousand-year period represented by our analysis.
Powder River has stored sediments in its active floodplain and three terraces. The uppermost and oldest terrace, the Kaycee Terrace, is typically 10-20 m above the modern riverbed. Leopold and Miller (1954) suggest an age range of 1600-4000 years before present (BP) for the Kaycee Terrace based on 14 C dates obtained outside the study area in the southern portion of the Powder River Basin within Wyoming. Below the Kaycee lies the Moorcroft Terrace, whose surface is 3.1-4.2 m above the modern riverbed, and the Lightning Terrace, the lowest and youngest of the Powder River terraces, which lies 2.1-3.4 m above the modern riverbed . The modern floodplain, the locus of the most active sedimentation along Powder River, is at most 1.76 m above the modern riverbed . Alluvial fan deposits are also present throughout the Powder River Basin. Alluvial fans typically extend only a short ( 1 km) distance from valley walls, with surfaces generally higher in elevation than those of the Kaycee Terrace.
Terraces identified along Powder River are not necessarily isolated from ongoing overbank deposition. At many locations along the river, the Lightning and Moorcroft terraces were fully inundated by flood waters in 1923 and 1978 , and likely also in 1890 (Metzger et al., 2020). During the 1978 flood, sediment deposited on the Lightning and Moorcroft terraces between Moorhead and Broadus represented 93% and 16% of the total suspended sediment transported during the flood at the Moorcroft stream gage  note that the suspended sediment brought into the reach at Moorhead was supplemented by erosion between Moorhead and Broadus, so the amount deposited can exceed the amount supplied by suspended sediment transport at Moorhead). No sediment was deposited on the Kaycee Terrace during the 1978 flood. Thus, on centennial timescales, overbank deposition remains an active process on the Lightning and Moorcroft terraces.
A wealth of geomorphic studies and supporting data provide an unusually extensive foundation for ongoing geomorphic research along Powder River. In 1975, Robert Meade established a series of cross-sections along Powder River that were surveyed at annual intervals in many of the ensuing years (Moody & Meade, 1990, 2020Moody et al., 2002). Cross-sections are labeled PR followed by the number of kilometers downstream from
This study also benefited from two additional sources of data. First, tree-ring data from 189 cottonwoods along 13 transects oriented along the axes of migrating meander bends were obtained by Schook et al. (2017). These data define the ages of modern floodplain deposits subject to erosion by Powder River. Second, in the Fall of 2016, a research-grade Light Detection and Ranging (LiDAR) survey of the study area was completed by the National Center for Airborne Laser Mapping (NCALM).

Overview
Given that the storage time distribution is equivalent to the age distribution of sediments leaving a storage reservoir (Bolin & Rodhe, 1973), the storage time distribution can be defined by dating sediments eroded during a specified time period, an approach initially proposed by Lancaster and Casebeer (2007) and later utilized by Lancaster et al. (2010). This method requires identifying a time period over which erosion is to be measured, quantifying the extent of erosion, and obtaining the age distribution of the eroded sediment using appropriate age-dating methods.
Powder River erodes sediment with ages ranging from less than a year to thousands of years, requiring multiple dating methods. Fortunately, at the beginning of this study, samples from Powder River's terraces had been dated using optically stimulated luminescence (OSL), and samples of the modern floodplain deposits were dated by Schook et al. (2017) using dendrochronology. The sampling strategy for these two dating programs was very different, and each requires a different approach for estimating the volume of eroded sediment. The results obtained using these disparate methods are integrated by analyzing geomorphic maps. The approach is summarized graphically in Figure 2.
To document locations and volumes of erosion, we selected the period from 1998 to 2013. This period was advantageous because cross channel topographic surveys were available, and because high-resolution aerial images were also available for these dates (additional details are provided in Schook et al. [2017]). The scale of the 1998 and 2013 aerial imagery is 1:40,000 and the resolution is approximately 1 m (Schook et al., 2017). Furthermore, channel positions for these years had already been digitized previously from aerial imagery (Schook et al., 2017), and these dates are also close in time to the 2016 aerial LiDAR survey, which provided detailed topographic data needed to identify different landforms along the river.

Measuring Meander Belt Width and Meander Wavelength
The width of the meadner belt can be related to valley confinement and also to the spatial and temporal scales of floodplain reworking. The meander belt is a region tangential to the outsides of meander bends (Chitale, 1970;Howett, 2017;Jefferson, 1902); this region was defined by drawing lines tangent to meanders in ArcGIS(v10.4.1). Meander belt width was estimated by drawing lines normal to the valley axis at 1 km intervals and averaging their lengths.
To assess confinement of Powder River's meanders by its valley margins, we also measured Powder River's meander wavelength (Camporeale et al., 2005). The average meander wavelength was estimated by first measuring the length of valley along its axis, and then dividing by two times the number of bends in our study reach (because a single meander wavelength includes two bends).

Geomorphic Mapping
Geomorphic mapping provides a foundation for integrating the results obtained during this study. Two maps were created, one documenting landforms of the valley of Powder River (Step 1 of Figure 2), and the other delimiting 10.1029/2021JF006313 6 of 21 sediments deposited by discrete episodes of lateral-accretion, which we term "lateral-accretion elements" (Step 2 of Figure 2). To identify different landforms, we relied on definitions and descriptions of Leopold and Miller (1954) and especially Moody and Meade (2008), who describe characteristic elevations of the modern floodplain and the Lightning, Moorcroft, and Kaycee Terraces. Terraces and alluvial fans were identified and defined using a variety of data sources, including USDA soil mapping, vegetation patterns observed on aerial imagery, and most importantly, aerial LiDAR (including first-return and ground-return data sets). The same data sources were also used to define the lateral accretion elements, which represent portions of the floodplain of approximately constant ages. The boundaries of lateral accretion elements are defined by topographic features such as floodplain scrolls, scarps bounding regions of differing elevations, vegetation patterns, and other discontinuities. Similar features have also been mapped by Durkin et al. (2018), who used them to define relative ages of preserved floodplain deposits.

Mapping Spatial Patterns of Erosion, 1998-2013
The active channel positions of Powder River in 1998 and 2013 were used to map spatial patterns of erosion. Schook et al. (2017) visually identified the active channel positions using aerial imagery (Step 3 of Figure 2), and these positions were used to determine the area of material eroded (Step 4 of Figure 2). The 1998 channel position was superimposed on the 2013 channel in ArcGIS. The portions of the 2013 channel exposed represent the area eroded between 1998 and 2013 ( Figure 3). The mapped lateral-accretion elements were extended to divide the areas of erosion into smaller erosion polygons, each representing eroded sediment of a nearly constant age (Step 5 of Figure 2; Figure 3).
The thickness of each erosion polygon is the difference between the elevations of its average upper surface and the average water surface (Step 6 of Figure 2). The elevation of the upper surface of an erosion polygon is defined as the average elevation of the adjacent bank as estimated from the digital elevation model (DEM) created from the 2016 aerial LiDAR data. Elevations were averaged within a 5-m (1/10 the channel width) wide band the length of the lateral-accretion element located five m away from the bank edge ( Figure 3). The five-m separation between the active channel and the sampling area was chosen to avoid very steep and irregular topography near the bank edge. This distance was chosen based on subjective field observations and assessments of the LiDAR data and aerial imagery. In the absence of detailed time series of bathymetric data, the lower surface of an erosion polygon is estimated by averaging the elevations of the water surface defined by the 2016 DEM in the area of each erosion polygon ( Figure 3). The volume of sediment eroded between 1998 and 2013 is obtained by multiplying the thickness of each erosion polygon by its surface area (Step 7 of Figure 2).
Each erosion polygon was assigned to either the active floodplain or to one of the three terraces based on the geomorphic mapping (Step 8 of Figure 2). The total volume of material eroded from each of these landforms was then computed (Step 9 of Figure 2), thus determining the contribution of each landform to the total erosion from 1998 to 2013.

Storage Time Distribution of the Modern Floodplain
The lowest surface within Powder River's valley is the modern floodplain. According to Moody and Meade (2008), the surface of the modern floodplain is "1.0-1.7 m above the riverbed, vegetated with cottonwoods (Populus deltoides, 0.01-0.1 m diameter) and with willows (Salix exigua, 0.01 m in diameter) interspersed with native and non-native grasses. Sedge grass (Carex sp.) forms a distinct vegetation band along the river edge near the break in slope between the bed and the bank and gradually disappears as the floodplain merges into a point bar." The modern floodplain is composed of material deposited within roughly the last 150 years making it suitable for cottonwood dendrochronology (Schook et al., 2017).
Plains cottonwood (Populus deltoides ssp. monilifera) trees typically germinate on newly deposited sediment on point bars on the insides of migrating meander bends along Powder River. The age of the cottonwood is thus a proxy for the age of the sediment on which it germinated (Everitt, 1968). As new floodplain sediments are deposited, the river migrates laterally, with younger cottonwood trees germinating on the newer sediment. Thus, tree ages document the progressive migration of the channel, providing a means of dating sediments as a function of distance along the axis of meander migration. The relation between age and the distance along the axes of Powder River's meander bends was determined by Schook et al. (2017). A smoothing spline was fit to these data, providing a means of estimating the ages of sediment within the interior of migrating meander bends along Powder River (Step 10 of Figure 2; Figure 4a) where dated trees were unavailable.
The axes of all meander bends within the study reach were drawn visually. After meander axes were drawn, ages were then estimated from the smoothing spline. By combining age estimates along meander axes with the mapped lateral-accretion elements, ages across most of the active floodplain could be defined (Steps 11 and 12 of Figure 2). To account for uncertainty in dating the active floodplain, data from individual erosion polygons were pooled into four age categories, each separated by 41.5 years (Step 13 of Figure 2). The size of these age categories (41.5 years) is the age of the oldest erosion polygon from the modern floodplain divided by four. Methods for quantifying uncertainty in ages and erosion volumes associated with each of these four age categories are discussed further below. This dating method is only valid if the direction of past meander migration can be inferred and where the history of migration has been unidirectional through time. In most meander bends, a distance is reached (as one progresses farther away from the active channel) where the geometry of the lateral-accretion deposits clearly records abrupt changes in the direction of meander migration, and at these transitions, the age of floodplain sediments is no longer a simple function of the distance along the meander axis. Fortunately, dates were obtained successfully for nearly all of the erosion polygons.

OSL Sample Collection and Dating
Samples were collected from nearly vertical banks over the 4-year period 2013-2016 at nine cross sections and at some additional sites. These were collected for OSL dating following standard sampling tube methods described by Nelson et al. (2015). We obtained OSL ages using standard methods designed for fluvial environments (Rittenour, 2008, and references therein). We extracted quartz from the samples following Nelson et al. (2015) and Gray et al. (2015) and adopted a Single Aliquot Regeneration protocol following Wintle (2000, 2003) with a preheat of 260°C. Aliquots were accepted with a rejection criteria of 10% for recycling ratio, 5% for recuperation, and 50% for equivalent dose and test dose error respectively. We performed dose recovery experiments and found our protocol can recover a given dose within 10%. We used these methods due to frequent problems with low quartz sensitivity and high recuperation rates that often led to low aliquot acceptance rates. Elemental concentrations were measured with in-house High Purity Germanium Gamma (HPGe) Spectrometry and Inductively Coupled Plasma Mass Spectrometry if the uncertainty on the HPGe data was larger than 50%. Dose rates were computed using the Dose Rate and Age Calculator (Durcan et al., 2015). To calculate ages, we applied the unlogged bootstrapped Minimum Age Model (Cunningham & Wallinga, 2012;Galbraith & Roberts, 2012) using functions in the Luminescence package for the programming language R (Kreutzer et al., 2020) and using an assumed value of 0.11 for σ b based on findings from Arnold and Roberts (2009) and Chamberlain et al. (2018).

Storage Time Distribution of Terrace Deposits
Age and erosion volume data for the Lightning and Moorcroft terraces were used to compute a storage time distribution for each landform. Volumes of erosion of the Lightning and Moorcroft terraces were calculated by determining the eroded area between cross-sections surveyed in 1998 and 2013. Eight cross-sections ( Figure 1a) with available OSL ages were included. Areas defined by the differences between surveys were treated as volumes by multiplying by a nominal width of 1 m in the downstream direction (Step 14 of Figure 2). The fraction contributed by the erosion of each cross-section to the total eroded volume of each terrace was computed, and an OSL age was assigned to each fraction (Steps 15 and 16 of Figure 2). If more than one OSL age was available from a particular eroding cutbank, then the ages were assigned to eroded sediments based on the elevation of each sample.
Cross-sections along Powder River do not adequately document erosion of the Kaycee terrace, because this landform is infrequently exposed along the river. As a result, samples from the Kaycee Terrace for OSL dating were not collected at cross-sections (Step 17 of Figure 2). Rather than distributing the ages based on direct observations of erosion (the procedure followed for all the younger landforms), Kaycee terrace ages were initially distributed based on rank (Step 18 of Figure 2) using the Weibull formula m/(n + 1) (Weibull, 1939), where m is the rank of a particular age (a rank of one is assigned to the youngest sample) and n is the number of samples. This approach assumes that the OSL dates of the Kaycee Terrace provide a complete and representative sample of the sediments eroded from this landform, which is unlikely to be true, and therefore we only use the median age of the Kaycee Terrace when we create the combined storage time distribution for all sediments eroded by Powder River from 1998 to 2013 (described below).

Combined Storage Time Distribution of All Eroded Sediments
The analyses described in Section 3.7 result in three separate data sets representing storage time distributions for modern floodplain, Lightning terrace, and Moorcroft terrace. As noted above, sediment eroded from the Kaycee terrace is represented by the median age of the available data (and its corresponding uncertainty). To create a single, combined storage time distribution for Powder River, the separate distributions (each covering a distinct and entirely separate age range) were scaled by the percentage each landform contributed to the total erosion from 1998 to 2013 (Step 19 of Figure 2) based on the areal landform mapping (Step 1 of Figure 2) and the analysis of spatial patterns of erosion related to lateral migration of Powder River (Steps 4-8 of Figure 2). The scaling process of combining the storage time distributions is illustrated graphically in the Results section.

Fitting Storage Time Distribution Functions
Exponential, Weibull, tempered Pareto, and the tempered first passage time distribution functions were each fit to the complete storage time distribution data using the curve fitting application within MATLAB (Version 2018a, Trust-Region algorithm; Table 1). The exponential distribution is a single parameter function that has been used to model well-mixed sediment storage reservoirs; it indicates that sediment within each age category has an equal chance of being eroded (Bolin & Rodhe, 1973;Everitt, 1968;Konrad, 2012). The other three distributions are included to assess the possibility that storage time distributions are heavy-tailed, with erosion probabilities weighted toward younger deposits in preference to older deposits. The Weibull distribution is a two-parameter function initially used to describe grain-size distributions (Rosin & Rammler, 1933) and subsequently applied in analyses of survival and failure rates (Lai et al., 2006). Moody (2017) successfully fit Weibull distributions to time-varying storage time distributions associated with a sediment pulse derived from hillslope erosion following a wildfire.
The tempered Pareto distribution has been shown to fit storage time distribution data from both modeling studies (Torres et al., 2017) and field data (Torres et al., 2020). Bradley and Tucker (2013) discuss the tempered first passage time distribution in a numerical modeling study of deposition and erosion by a meandering river system. The Pareto and first passage time distributions are "tempered" when some process such as a physical boundary or other limit constrains the maximum possible age of stored sediment. For both distributions, "tempering" is achieved by adding an exponential term scaled by the maximum possible age (Table 1).
In addition to the statistics provided directly by the MATLAB curve-fitting application (e.g., adjusted coefficient of determination and root mean square error), we also sought an additional means of assessing the goodness-offit of these functions to our data. Remarkably, it turns out that the three heavy-tailed distributions are strongly tempered when fit to our data, and are all essentially equivalent to the exponential distribution (note that the "tempering" maximum age scale is determined by the curve-fitting routines rather than being imposed from the data before curve fitting). This interpretation (fully documented below in Section 4) allows us to limit our assessment of goodness-of-fit to the exponential distribution.
We used the Anderson-Darling test (Feigelson & Babu, 2021;Stephens, 1986) implemented in the MATLAB function adtest, to assess the goodness-of-fit of the exponential distribution to our observed storage time distribution. We used the Monte-Carlo procedures in adtest, with a tolerance of 0.01 to compute p-values.

Quantifying Uncertainty
The storage time distribution is based on estimates of eroded sediment volumes and their ages. We assessed uncertainty in as many of our measurements as possible, and propagated error estimates through all our computations (Topping, 1972).

Uncertainty in Locating Boundaries From Aerial Imagery
Polygons of lateral erosion were identified from the 1998 and 2013 mapping of river channel boundaries. These polygons are used to assign eroded sediment to different landforms (Step 8 of Figure 2), and also to assign eroded sediment to discrete time periods in the history of the active floodplain using the mapping of lateral-accretion elements (Step 5 of Figure 2). To evaluate the uncertainty in these estimates the locations of static features along the river were mapped on the aerial imagery at different times. For example, surveyed cross-sections demonstrate that the left bank (facing downstream) of Powder River at the location of the USGS Moorhead stream gaging station has been stable for decades (Moody & Meade, 2018), so this site was repeatedly mapped on aerial imagery. Uncertainty in locating this boundary at different times on aerial imagery averaged about ±5 m, so this value was assumed to represent uncertainty associated with any estimate of boundary position as mapped throughout the study area on the 1998 and 2013 aerial imagery.

Uncertainty Estimates for the Thickness of Eroded Polygons
The thickness of each erosion polygon is defined as the difference between upper and lower surfaces. The average elevations of both surfaces were estimated by averaging data from the DEM constructed from the NCALM Li-DAR data. To evaluate the uncertainty associated with each estimate of average elevation, the standard deviation was computed. Standard deviations for the upper surface of the eroded polygons ranged from 0.03 to 2.01 m, while standard deviations for the lower surface ranged from zero to 0.76 m.

Uncertainty of Binned Modern Floodplain Storage Ages
The age of each erosion polygon within the modern floodplain was determined by the spline function fit to the ages of cottonwood trees and their distances along meander axes (Figure 4) Note. Columns for λ and k provide best-fit values from the MATLAB curve-fitting routines, with 95% confidence intervals (C.I.) in parentheses. The column "Number within error bounds" presents the number of points for which the best-fit equation passes within the uncertainty envelope of the data. The p-value is obtained from the Anderson-Darling test implemented by the MATLAB routine adtest.

Table 1 Cumulative Distribution Functions Fit to Powder River Storage Time Distribution Data
standard error of ±18.54 years; this value was used to represent the age uncertainty of all polygons on the active floodplain. The average age uncertainty for each age bin is a result of summing the age uncertainties in quadrature. Age uncertainty associated with the four age bins ranged from 24 to 82 years.
The binning process creates an additional source of uncertainty related to the possibility that erosion polygons could be assigned to the wrong age bin. Several steps were followed to quantify this uncertainty. First, the distance along the x-axis was divided into five segments defined by the intersections between the spline function and the boundaries of the age categories (Figure 4). Within each of the five distance segments, the number of trees in each of the four age categories was counted. The uncertainty associated with assigning erosion polygons to a given age category (as determined by the spline curve) was then defined as the fraction of the number of trees in all the other (unassigned) categories combined. For example, the distance increment from 112 to 349 m is assigned to the age category 41.5-83 years ( Figure 4). Within this distance increment, there are 52 trees within this age range and 17 trees with ages outside this age range, 10 younger than 41.5 years and 7 older than 83 years. This leads to an uncertainty estimate of 17/69 or 24.6% for assigning sediment to this age category within this distance increment (Figure 4). Uncertainties range from 21% to 44% (Figure 4). These uncertainty estimates were added to the age uncertainty discussed above and then propagated through all computations.

Uncertainty Estimates for Terrace Volumes and Ages
In computing the storage time distribution of the eroded terrace deposits, uncertainty estimates are needed for ages and volumes of erosion. Age uncertainties are directly provided by the uncertainties associated with OSL age dating (Table 2). Uncertainties associated with mapping erosion volumes of the terraces from aerial imagery have already been discussed. These age and volume uncertainties were propagated through all computations.

In
Step 14 of Figure 2, the age distribution of erosion volumes from the Lightning and Moorcroft Terraces was assessed from eight surveyed cross-sections (in addition to the assessment of the total volume of erosion from each terrace using historical aerial imagery, i.e., Step 9 of Figure 2). The survey data themselves are very precise, with elevation uncertainties of 0.01 cm and horizontal uncertainty of 0.1 m (Moody & Meade, 2018), and therefore the surveyed erosion volumes appear to have insignificant uncertainty compared to other sources of uncertainty, so uncertainty associated with surveying has been neglected.
Additional uncertainty can be associated with the use of a limited number of cross-sections to determine the age distribution of terrace sediments eroded between 1998 and 2013. Ideally, age data from each terrace would be distributed randomly throughout the eroded sediment. However, a random sampling protocol was impractical given the resources available to us, and as a result, errors associated with non-random sampling of terrace deposits cannot be rigorously assessed.

Uncertainty Associated With Mapping Terraces
Geomorphic mapping (Step 1 of Figure 2) provides the basis for allocating eroded volumes to each of the four alluvial landforms of Powder River (Step 19 of Figure 2). To assess uncertainty created by errors in geomorphic mapping, we compared maps created during this study with similar maps of Moody and Meade (2008, their Figure 2) showing the locations of the modern floodplain and the Lightning, Moorcroft, and Kaycee Terraces within a subset of our study area. We overlaid the erosion polygons for the area covered by both sets of maps and determined the total areas of erosion associated with each landform based on each set of geomorphic maps, resulting in two entirely independent estimates of the fraction of the total erosion associated with each landform. We computed the absolute value of the differences between estimates of the fractional contribution to the total erosion for each landform obtained from the two different geomorphic maps. These differences represent the uncertainty associated with geomorphic mapping for estimating the erosion contributed by each landform to the total erosion between 1998 and 2013. Uncertainties associated with geomorphic mapping were added to other sources of error and propagated through all computations.   Camporeale et al. (2005)

for Unconfined Meanders
Our method for assessing uncertainty in geomorphic mapping is particularly rigorous because the two maps were created using different data covering different time periods. In this study, geomorphic maps were based on high resolution aerial imagery, aerial LiDAR, soil maps, and other remotely sensed data available through 2017, while the maps created by Moody and Meade (2008) are based on surveyed cross-sections, field observations, and aerial imagery largely representing conditions immediately after 1978 flood. Because of differences in resolution and time period, the comparison of our mapping with that of Moody and Meade (2008) could overestimate uncertainties associated with our landform mapping.

Meander Morphology and Geomorphic Mapping
A total of 47 bends were counted in the study reach over a downvalley distance of 17 km. Thus, the meander wavelength of Powder River is 720 m. The meander belt of Powder River averages 900 m wide, with a standard deviation of 150 m (Figure 5a). Camporeale et al. (2005) found that the ratio of meander belt width to meander wavelength can indicate whether a river channel is confined or unconfined, with values 3.4 characteristic of confined rivers, and values 3.4 unconfined. Here the ratio of meander belt width to meander wavelength is 1.25; thus, the study reach of Powder River represents confined meandering.
The modern floodplain is generally adjacent to the active channel except where terraces project out from the valley margin (Figure 5a). The Lightning, Moorcroft, and Kaycee terraces tend to be progressively farther from the active channel, and typically only meet the active channel on the outsides of meander bends . The Kaycee terrace and alluvial fan deposits are primarily preserved along the valley margins. The modern floodplain and Lightning and Moorcroft terraces are generally contained within Powder River's meander belt, while the Kaycee terrace and alluvial fans generally lie outside the meander belt (Table 2). This suggests that the modern floodplain, Lightning, and Moorcroft terraces are all deposits that are "frequently" reworked by meander migration of Powder River. If the meander belt is extended to 3.4 times the meander wavelength, as would be expected for unconfined meandering (Camporeale et al., 2005), then this unconfined meander belt would include nearly all mapped landforms (Table 2), encompassing the Kaycee terrace, alluvial fans, as well as portions of the valley margins.
The mapped lateral-accretion elements average approximately 40-m wide, are typically lenticular, and match the curvature of the channel at the time of deposition ( Figure 5b). As the distance from the channel increases the shapes become more irregular. Shifting channel migration patterns through time are revealed by cross-cutting relations between boundaries of lateral-accretion elements (Figures 3 and 5b).
The majority of erosion from 1998 to 2013 was on the outsides of meander bends. Erosion polygons typically form thin slices whose widths are a fraction of the channel width (red areas in Figure 5c).
Between 1998 and 2013, nearly 7.0 × 10 5 m 3 of sediment was eroded from the approximately 17 km study reach. The modern floodplain contributed more sediment (39% ± 9%) than any of the terraces individually, but the three terraces in total contributed 56% more sediment (61% ± 9%) than the modern floodplain. Contributions from individual terraces decreased with increasing age: erosion from the Lightning, Moorcroft, and Kaycee terraces represented 28% ± 7%, 20% ± 5%, and 13% ± 3%, respectively. Alluvial fan deposits provided only 1.5% of the material eroded over the 15-year period. Due to the minimal input of material eroded from alluvial fans and lack of age constraints for these deposits, input from alluvial fans was not included in the assessment of the storage time distribution.

Storage Time Distribution
The storage time distribution includes dates from four landforms: the modern floodplain and the Lightning, Moorcroft, and Kaycee terraces. Ages of erosion polygons for the modern floodplain range from 24 to 144 years ( Table 3). The 186 erosion polygons contained within the study reach represent a total of 1.47 × 10 5 m 3 of eroded sediment or approximately 55% of the sediment eroded from the entire modern floodplain. Age uncertainties ranged from ±24 years to ±82 years, while volume uncertainties are between 2% and 12%. Age uncertainties are primarily associated with the use of the spline function to place erosion polygons into the appropriate age category. The Lightning terrace ages range from 400 to 950 years with 2 sigma uncertainties ranging from ±70 to ±190 years (Figure 6a). The Moorcroft terrace ages range from 1130 to 1897 years with 2 sigma uncertainties ranging from ±96 to ±156 years. Ages for the Kaycee terrace range from 1610 to 6100 years with 2 sigma uncertainties ranging from ±330 to ±1370 years. The median age of the Kaycee Terrace is 5000 ± 680 years.
When the four distributions of Figure 6a are scaled according to the erosion from each landform and combined, a single storage time distribution is obtained for the entire study area that extends from the present to 5,000 years ago (Figure 6b). The mean of the storage time distribution (the residence time) is 928 ± 167 years, while the median is 774 ± 111 years. About 70% of the eroded sediment is younger than 1000 years. The combined distribution of Figure 6a shows a discontinuity as an abrupt horizontal offset between the modern floodplain and the Lightning Terrace ages of approximately 250 years (Table 3; this feature of our data is addressed in Section 5).     (Table 2).

Storage Time Distribution Functions
The four mathematical functions fit to the observed storage time distribution appear to be almost identical (Figure 7). Coefficients of determination (r 2 ) are generally high, ranging from 0.76 for the tempered first passage time distribution to 0.87 for the tempered Pareto distribution (Table 1). Root-mean-square errors range from 0.09 (tempered Pareto distribution) to 0.13 (tempered first passage time distribution). The Anderson-Darling goodness-of-fit statistic for the exponential distribution of 0.003 is significant (Table 1).
The similarity between all four fitted storage time distributions functions suggest that the data are exponentially distributed, and that the heavy-tails of the Weibull, tempered Pareto, and tempered first passage time distribution are poorly defined. This can be readily demonstrated by examining fitted parameter ranges for the non-exponential terms in these functions. For the Weibull distribution, the 95% confidence interval for the fitted value of the exponent k includes 1; for this value of k, the Weibull distribution reduces precisely to the exponential distribution. For the tempered Pareto distribution, the 95% confidence interval for k includes the value of 0, which eliminates the power term, and the tempered Pareto distribution also becomes equivalent to an exponential distribution. Results for the tempered first passage time are somewhat more complex, but lead to a similar conclusion. For small values of x, the argument of the error function in the tempered first passage time distribution is large, leading a value of 1 for the error function, and the tempered first passage time distribution becomes an exponential distribution. For large values of x, the argument of the error function is small, and the error function becomes important, but now the exponential term is close to 0, and the resulting value of the function is close to 1, as would be given by the exponential distribution. Thus, all the functions are essentially exponential. Furthermore, the exponential scaling factor λ values are very similar for all the fitted functions, ranging from 785 to 1274 years. Statistically, then, all of these fitted equations are exponential, with similar exponential scaling factors.

Discussion
The approach we have described determines the complete distribution of sediment storage times for Powder River. Our method quantifies erosion volumetrically over a given time period, as well as the depositional age of the eroded sediment. We believe that we have described the only field method available for estimating a complete storage time distribution for alluvial rivers. The results confirm previous research (Bradley & Tucker, 2013;Meade, 2007;Pizzuto, 2020) suggesting that sediments deposited on floodplains may remain in storage for centuries to millennia ( Figure 6). Sediments stored on the modern floodplain of the contemporary Powder River were stored for less than 166 years and comprise 39% of the eroded sediment. The remaining 61% of eroded sediments had been stored in terraces for up to 5,000 years.

Confinement of Powder River's Meander Belt by the Kaycee Terrace and Alluvial Fans
The meander belt represents the area of Powder River's valley actively reworked by meander migration. On Powder River, the meander belt width is 1.25 times the meander wavelength, considerably narrower than the value of 3.4 meander wavelengths cited by Camporeale et al. (2005) as characteristic of unconfined meandering rivers. This suggests that meandering of Powder River is laterally confined.
Landforms primarily confined within Powder River's meander belt include the modern floodplain, Lightning, and Moorcroft terraces, while the Kaycee terrace and alluvial fans are mostly located outside the meander belt (Table 2) and represent confining elements of the meander belt along Powder River's valley margin (Fryirs et al., 2016). The river's modern floodplain, Lightning, and Moorcroft terraces all increase their elevation through vertical accretion on centennial timescales , while the Kaycee is not subject to episodic flooding and aggradation. These observations all suggest that within the timescale of formation of the meander belt, the modern floodplain, Lightning, and Moorcroft terraces represent active alluvial sediment reservoirs that continue to exchange sediment with Powder River. Whereas, sediments stored in most of the Kaycee Terrace segments (Figure 5a) represent alluvium and colluvium stored during an earlier stage in Powder River's evolution.

Scaling and Sediment Budget Estimates of Mean Storage Timescales
Results of Figure 7 and the parameters fitted to the four storage time distributions all indicate that the storage time distribution is exponential, representing a sediment storage reservoir defined by the river's meander belt that is "well-mixed," such that all deposits are equally likely to be remobilized by Powder River's lateral migration. Furthermore, best-fit and 95% confidence intervals for the exponential "timescale" parameter λ, with values of 785-1274 years and a range of 230-2300 years (Table 1), can be related to the timescale required for Powder River to rework its meander belt. This value is approximately 900 years given the meander belt width of 900 m divided by Powder River's approximate meander migration rate of 1 m/yr. These results provide a mechanistic explanation for the observed storage time distribution by relating the exponential mean timescale to the physical process of sediment reworking by lateral migration.
A second independent assessment of Powder River's mean storage timescale can be obtained by analyzing erosion rates and storage volumes, though additional assumptions are required. For a steady-state system in which rates of erosion and deposition are equal and constant with time, the residence time (mean storage time) is given by the ratio of the stored volume to the annual volumetric erosion (or deposition) rate (Bolin & Rodhe, 1973).
Here we assess this ratio using areal rather than volumetric estimates, because the vertical extent of erosion is more uncertain than areal measurements. The total areal erosion from 1998 to 2013 illustrated in Figure 4 is 2.9 × 10 5 m 2 , while the total area within Powder River's meander belt is the product of the meander belt width (900 m) and the valley length of our study reach (17 km), equivalent to 1.5 × 10 7 m 2 . Dividing the total area by the annual reworking rate (approximately 1.9 × 10 4 m 2 /yr) gives a storage timescale of 780 years, a value that agrees remarkably well with the range of best-fit values obtained from the storage time distribution of 785-1274 years and that is also consistent with the estimate of 900 years obtained above from the ratio of meander belt width to lateral migration rate.
These results illustrate how complete storage time distributions for confined meandering systems can be obtained solely from simple measurements from aerial imagery. However, these methods are strictly valid when fluvial systems have maintained constant and approximately equal rates of erosion and deposition through the period represented by the stored sediment (so that recent erosion rates can be extrapolated throughout the entire time required to rework the meander belt).

Implications of an Exponential Storage Time Distribution
An exponential storage time distribution for deposits of confined meandering rivers has important implications for routing of sediments, associated contaminants, and particulate organic carbon. If storage time distributions are exponential, the mean storage time (residence time) provides a useful metric for assessing storage timescales, whereas if storage time distributions are heavy tailed, the mean may not exist, and a small percentage of stored sediment may persist for exceptionally long timescales.
An exponential storage time distribution also allows for simple computations of sediment travel times through a reach. This arises when the sediment budget is known and can be treated as approximately constant through time and throughout the reach of interest. When these conditions are met, the complete travel time distribution for sediment and associated conservative contaminants through the reach can be computed analytically, without requiring a complex numerical model. The approach is summarized by Pizzuto et al. (2022) and Torres et al. (2017), but cannot be implemented for Powder River because an accurate sediment budget is not yet available. Nonetheless, the existence of an analytical method for computing sediment travel times provides an important means of readily assessing the influence of storage on sediment delivery.

Additional Sources of Uncertainty
Significant uncertainties remain in our estimates of deposit ages and fractional volumes. These uncertainties often approach ±10% ( Figure 6, Table 3). As a result, our methods may fail to accurately detect small amounts of older sediments eroded from the Kaycee Terrace, and these deposits, if accurately sampled and dated, could conceivably support the interpretation of a heavy-tailed storage time distribution for Powder River. However, given the precision of available dating methods and the uncertainties associated with geomorphic mapping and sampling, it may not be possible to identify long-tails in the storage time distribution that arise from erosion of a few older sediments.
It is also possible that some of our results have been biased by our sampling methods. A longer study reach might have allowed for sampling a greater number of exposures of the Kaycee Terrace, which could potentially have resulted in larger volume of sediment that could have defined a heavy-tailed storage time distribution. The use of the water surface instead of the streambed as the lower boundary for erosion could also be considered a limitation of our study, as our approach cannot account for sediment eroded from areas submerged by the river. However, while our approach likely underestimates the total erosion, the distribution of ages of the eroded sediment, the key variable for our interpretation, is probably unaffected by our use of the water surface as a datum. Furthermore, detailed time series of bathymetric data for Powder River that could be used to define the full vertical extent of erosion do not exist.
The timescale associated with the exponential storage time distribution appears to be well-correlated with Powder River's meander belt width and the current meander migration rate of approximately 1 m/yr. However, meander migration rates have varied systematically since 1830, decreasing significantly over the last 40 years, changes that have also been accompanied by an increase in Powder River's sinuosity and a decrease in its width (Schook et al., 2017). By interpreting the storage timescale in terms of meander migration rates, we have extrapolated the average migration rate for 1830-present to cover the entire timescale of our age dating, encompassing at least 5,000 years. Meander migration rates have likely varied systematically over this time period due to changes in climate and other drivers. These effects cannot be evaluated rigorously at present, but could represent a source of error.
It is also interesting to speculate on how storage distributions might vary for different grain size fractions transported by Powder River. Stored sediment in meandering rivers is typically divided into two categories, sand-and gravel-sized point bar deposits, mostly deposited from bedload, and finer overbank sediments deposited from suspension (Ghinassi et al., 2018). While no field studies are available to assess differences in storage time distributions associated with these different grain size fractions, a modeling study by Ackerman and Pizzuto (2016) suggests that the storage time distributions of overbank and point bar deposits are similar. This interpretation likely applies to Powder River because coeval deposits represented by the lateral accretion elements mapped in Figure 4 include both point bar and overbank deposits, and these sediments would also both be remobilized at the same time when Powder River's channel migrates laterally and removes them. This hypothesis, of course, should be tested by additional studies.

Conclusion
Powder River's meander belt is considerably narrower than the 3.4 meander wavelengths associated with unconfined meandering, suggesting that Powder River's meanders erode and deposit sediment within a laterally confined valley. The modern floodplain, Lightning terrace, and Moorcroft terrace are all primarily located within Powder River's meander belt, and thus these landforms represent the active storage reservoirs exchanging sediment with Powder River. Landforms confining active meander migration include alluvial fans and the Kaycee terrace.
Ages of sediment eroded between 1998 and 2013 from a 17 km length of Powder River's valley are exponentially distributed, with a best-fit mean age of 824 years (95% C.I. 610-1030 years). We also obtained two independent estimates of the mean age of stored sediment to support our analysis of the age distribution of eroded sediment. The ratio of meander belt width (900 m) to time-averaged meander migration rate (∼1 m/yr) is 900 years, while the ratio of meander belt area to the annual areal erosion is 780 years. Both of these independent estimates of the mean storage time agree remarkably well with the mean value of the exponential distribution extracted from our data.
Our finding that Powder River's storage time distribution is exponential does not support our initial hypothesis that floodplain storage time distributions are heavy-tailed. We speculate that heavy-tailed distributions are more likely to be found where meandering is unconfined.
If the results of this paper can be generalized, analyses of meander belt storage reservoirs will be greatly facilitated. The scaling parameter of the exponential distribution describing both the storage time and age distributions can be readily estimated from either of two ratios: the meander belt width divided by the meander migration rate, or alternatively, the meander belt area divided by the annual areal erosion rate. These ratios can be estimated from aerial imagery without the extensive field sampling, sediment dating, and landform classification methods used here. Defining the storage time distribution can be helpful to better understand sediment travel times through a reach, and if a sediment budget is available, relatively straightforward analytical methods are available to compute the complete travel time distribution, results that are invaluable for estimating the downstream migration of sediment-associated contaminants, the fate of particulate organic carbon, the timescale required to realize benefits from sediment-related restoration strategies, and the delivery of climate and tectonic driven sedimentary signals to sedimentary basins.