A Bayesian mixing model framework for quantifying temporal variation in source of sediment to lakes across broad hydrological gradients of floodplains

Paleolimnological reconstructions provide insights into hydrological variability of dynamic floodplain lakes. However, spatial and temporal integration of multiple reconstructions often remains underdeveloped because the efficacy of different paleolimnological measurements varies among lakes due to gradients in energy of floodwaters and sediment composition. Here, we use linear discriminant analysis to identify 10 significant elemental concentrations in sediment obtained from multiple sampling campaigns that distinguish among three end‐member allochthonous sources for lakes in the Peace‐Athabasca Delta (PAD; Canada): Athabasca River, Peace River, and local catchment runoff. Over 90% of sediment samples were correctly classified into original groups after cross‐validation due to the distinctiveness of the three end‐members, which permitted development of a robust Bayesian mixing model to discern the relative contributions of sediment from the three sources. We evaluate performance of the mixing model via application to sediment cores from two adjacent lakes in the Athabasca sector of the PAD and demonstrate its effectiveness to discriminate three known hydrological phases during the past 300 years. Notably, model results indicated that ~ 60% of the sediment originated from the Peace River during the largest ice‐jam flood event on record (1974), which was unrecognized by other methods. The approach provides a new, universal method that can be applied across the full range of sediment composition to quantify changes in source, frequency, and magnitude of sediment delivery by river floodwaters to lakes and is transferable to other dynamic floodplain landscapes where broad range of sediment composition challenges application of other approaches.

locations receive less frequent and weaker influence of floodwaters, and their sediments are typically organic-rich. Consequently, analysis of biological (e.g., diatoms, plant macrofossils) and organic geochemistry (e.g., carbon and nitrogen elemental and isotope composition) components are more effective to reconstruct past hydroecological change and variability (Wolfe et al. 2005(Wolfe et al. , 2008bBrock et al. 2010;Kuerten et al. 2013). Integration of paleolimnological records from multiple lakes that span the broad gradients of frequency and intensity of flooding may provide important additional insight into factors regulating hydroecological conditions across floodplain landscapes that are unavailable from analysis of one or a few lakes at either end of these energy gradients. However, reconciling records based on physical methods at some sites and biological and geochemical methods at other sites presents challenges, which frequently constrain reconstructions to qualitative or semiquantitative descriptions of past hydrological conditions (Wolfe et al. 2008a;Velez et al. 2012;Kuerten et al. 2013). Furthermore, analyses of most physical, biological, and geochemical variables in sediment cores do not inform about the source of floodwaters, which may be important to discriminate when more than one river has potential to be the source of flooding. Thus, development of a method that allows consistent measurement of the same variable in cores from lakes across the gradient of mineral-rich to organic-rich sediment and is capable of fingerprinting the source of sediment would represent a useful advance for elucidating drivers of hydroecological change in floodplains.
Fluvial sediment transported to floodplains possesses materials derived from potentially different regions within watersheds due to variation of the surficial geology, and these mixtures are deposited into lakes during flood events where they accumulate vertically and form a stratigraphic profile (Kitch et al. 2019). Sediments deposited into floodplain lakes, thus, provide an integrated temporal record of the parental sediment sources, each with their own geochemical fingerprint (Collins et al. 2010(Collins et al. , 2017Blake et al. 2018). Sediment fingerprinting by trace element concentrations or isotope ratios has been used to discriminate among parental sourcematerials (Walling and Woodward 1995;Hughes et al. 2009;Lamba et al. 2015a,b;Stewart et al. 2015). Unlike other measurements that vary in their efficacy across gradients of mineral and organic content, trace element concentrations can be quantified reproducibly across the spectrum from mineral-rich to organic-rich lake sediment, as occurs along gradients of river connectivity within floodplains and deltas (Lintern et al. 2016;Kay et al. 2020;Klemt et al. 2020). Thus, elemental concentrations may be ideal input tracers to track changes in the sources of sediment supplied to lakes across hydrologically complex floodplains, deltas, and other lake-rich landscapes where it is important to discriminate past variation of multiple sources of sediment, including floodwaters provided by one or more possible rivers. Indeed, mixing models have been used to quantify the role of anthropogenic disturbances within upland portions of watersheds on sediment deposition to downstream lakes (Cooper et al. 2015;Lizaga et al. 2019;Wynants et al. 2020). For instance, Wynants et al. (2020) recently used a Bayesian mixing model (BMM) framework and measurements of trace elements in lake and river sediment samples at strategically selected locations within Lake Manyara (Tanzania) and its many tributaries to quantify geochemical differences of river sediment supplied to the lake, and to reconstruct temporal variation in sediment sources from stratigraphic profiles. The approach revealed that a recent increase in siltation was linked to expansion of agricultural activities within the watershed of one of the lake's many tributaries-knowledge that assisted land management planning and protection of aquatic habitat and water quality of this large East-African Rift Lake. Defining (using concentrations of elements) and quantifying (using a BMM) the relative importance of sediment sources, and the temporal variability of these sediment sources, to lakes in floodplains, deltas and lake-rich landscapes has the potential to enhance quantitative paleohydrological reconstructions, and better inform water resource management.
The approach elegantly demonstrated by Wynants et al. (2020) presents new opportunity to further probe the hydrological evolution of the Peace-Athabasca Delta (PAD) in northeastern Alberta (Canada), a site of intensive paleolimnological investigation for the past 20 yr Kay et al. 2019Kay et al. , 2020. The PAD is a 6000 km 2 boreal freshwater floodplain where episodic input of floodwaters caused by ice-jams on the Peace and Athabasca rivers recharges its abundant shallow, perched lakes. Several upstream and resistance forces may contribute to the occurrence of an ice-jam flood event Conly 1998, Beltaos andBonsal 2021;Lamontagne et al. in  Paleolimnological studies have demonstrated the decline in magnitude and frequency of ice-jam flood events and lakelevel drawdown began in the early 20 th century and have provided key insight into the paleohydrological record of individual lakes. However, spatial integration is not yet fully developed because different paleolimnological measurements have been used among lakes depending on the strongly varying mineral-rich to organic-rich sediment composition across the strong hydrological gradients of the delta.
Recent analysis of elemental concentrations in sediments of lakes across the PAD has demonstrated significantly higher concentrations of some elements in sediment conveyed by the Peace River compared to sediment sourced from the Athabasca River, which is attributed to the different underlying geological terrane of their catchments (Owca et al. 2020). This finding affords opportunity to develop a novel approach capable of quantitative reconstruction of sediment supplied by the Peace River vs. Athabasca River to lakes, which is important in the vast low-relief PAD where multiple factors have potential to differentially alter flood regimes of these rivers, and large central areas of the delta may variously receive inflow from both rivers. In this study, we capitalize on differences in elemental concentrations in sediment supplied by the Peace River and Athabasca River, and from the local catchment in the absence of river floodwater, to evaluate the utility of a BMM (MixSIAR; Stock and Semmens 2016;Stock et al. 2018) in quantifying relative source contributions from these three end-members. We apply a MixSIAR model to stratigraphic profiles of elemental concentrations from sediment cores of two lakes in the PAD, where knowledge of past marked hydrological change and sediment composition has been previously obtained using more conventional methods (Wolfe et al. 2008b;Kay et al. 2019), to assess performance of the model. As we demonstrate, the methodology fosters opportunity for improved characterization and quantification of hydrological change over time and space in complex floodplain environments, such as the PAD, and is readily transferable to other comparable water-rich landscapes. This knowledge is important for stewardship decisions and mitigative actions to maintain abundance of water and ecological integrity.

Materials and procedures
Three steps were employed to develop, apply, and evaluate the methodological approach. First, samples of surficial lake sediment, river bottom sediment, and riverbank sediment were obtained from an array of sites that span the hydrological gradients of the PAD, and each sample was classified into one of three distinct a priori groups of allochthonous sediment input (hereafter referred to as end-members): (1) Athabasca River sediment, (2) Peace River sediment, and (3) local catchment. These three end-members represent the three major pathways of sediment delivery to lakes in the PAD under flood (Athabasca River, Peace River) and nonflood (local catchment) conditions and are likely to be preserved in stratigraphic records. Samples were designated to an end-member based mainly on their geographic location, lake water isotope composition, and lake water chemistry data. Second, linear discriminant analysis (LDA) was conducted on sediment samples to select a subset of elements that best discriminated among the three end-members (via removal of collinear elements) and to determine if the elemental composition of each end-member differed significantly. Cross-validation was then performed to evaluate ability of the LDA to discriminate endmember sediment sources in the PAD based on differences in elemental composition of lake sediment. Third, means and standard deviations (SDs) of the elements, which differed significantly among Peace River, Athabasca River, and local catchment sediments (as identified by the LDA), were entered into a MixSIAR model to represent each of the three end-members. The MixSIAR model was then applied to stratigraphic profiles of element concentrations obtained from two lakes (PAD 30 and PAD 31) known to have experienced three distinct hydrological phases during the past 300 years and span a broad range of sediment composition (40-90% mineral matter; Wolfe et al. 2008b;Kay et al. 2019). This last step provided opportunity to evaluate performance of the model with independent data. Further details for each step are provided below.

Collection, analysis, and designation of sediment samples to the three end-members
During 2017-2019, samples of lake surface sediment, surficial river-bottom sediment, and riverbank sediment were collected in the PAD (Fig. 1). Surface sediment samples were obtained using a mini-Glew gravity corer (Glew 1991) from a network of 44 lakes across the PAD in September 2017 (8 flooded by the Athabasca River, 3 flooded by the Peace River, 33 local catchment), and from a subset of 13 of those lakes in summer 2018 after they received floodwaters from an ice-jam flood event during that spring (6 flooded by the Athabasca River, 7 flooded by the Peace River; Owca et al. 2020). Nine river-bottom samples were obtained in 2019 from the Peace River using a Petite Ponar grab sampler. No river-bottom samples were obtained from the Athabasca River. Riverbank sediment samples were collected from exposures along the Peace River (n = 10) and Athabasca River (n = 36) during summer 2019. The samples of lake surface sediment, river-bottom sediment, and riverbank sediment were homogenized, and subsamples were analyzed for a suite of elements at ALS Canada (Waterloo) following EPA method 202.2/6020A, which uses an aqua regia digestion.
To develop a geochemical fingerprint for each major source pathway, the sediment samples were assigned to one of the three end-members (Athabasca River sediment, Peace River sediment, or local catchment). Riverbank and river-bottom sediment samples were assigned to end-members of the river from which they were collected. Flooding in the PAD is variable between years, and lakes can be flooded by at least one of the Peace River or Athabasca River depending on the lake's geographic location and elevation (Remmer et al. 2020a). Lake surficial sediment samples collected in 2017 were assigned to the Peace River and Athabasca River end-members based on the sampling location, and determination of whether the sample likely incorporated inputs from river floodwaters was based on evidence from repeated seasonal measurements of lake water isotope composition and water chemistry during 2015-2019, following Remmer et al. (2020a). Surficial sediment samples from the nonflooded lakes were placed in the "local catchment" end-member. This end-member captures lake sediment signatures in the absence of recent flooding and includes fluvio-deltaic sediment eroded from the catchment. The 1-cm thick surficial sediment samples obtained from the lakes in 2017 likely captured sediment deposited within the prior 2-5 yr based on our experience dating cores from > 30 lakes in the PAD Kay et al. 2019). Lakes located in the central region of the delta were excluded from the study because sediment could be derived from multiple sources, which potentially confounds attempts to characterize the geochemical fingerprint of each of the three individual sources. Lake surficial sediment samples collected in 2018 were exclusively from lakes flooded in spring of that year and were assigned to the Peace River and Athabasca River end-members based on geographic location, following Remmer et al. (2020b) and Owca et al. (2020). Only the uppermost distinctive gray flood layer was collected to ensure that samples obtained exclusively represented sediment deposited during the spring ice-jam flood event. Lakes that were sampled in both 2017 and 2018 were allocated to an end-member in each year based on the above criteria and treated as independent samples in subsequent analyses.
In summary, we allocated sediment samples to the three end-members as follows. (1) The Athabasca River end-member is represented by the surface sediment samples from lakes in the southern Athabasca sector of the PAD considered to have recently received floodwater from the Athabasca River, plus the Athabasca River riverbank samples (total for end-member = 50).
(2) The Peace River end-member is represented by the surface sediment samples from lakes in the northern Peace sector of the PAD considered to have recently received floodwater from the Peace River, plus the Peace River surficial river-bottom and riverbank sediment samples (total for endmember = 29). (3) The local catchment end-member is represented by surficial sediment samples from lakes in both the Peace and Athabasca sectors that were deemed to have not recently received river floodwaters (total for endmember = 33).

Linear discriminant analysis
Key requirements for successful application of mixing models include discrimination among end-members, and that input tracers are spatially and temporally conservative in their environmental behavior (Motha et al. 2002;Blake et al. 2018). This implies no changes in the composition of the tracers within sediment profiles before or after deposition (Koiter et al. 2013;Belmont et al. 2014;Laceby et al. 2017). Thus, prior to performing LDA, we removed elements from the initial suite known to exhibit nonconservative behavior (i.e., elements sensitive to redox reactions, mobile in sediments, or readily altered by diagenetic processes) from the analysis, as recommended by Smith et al. (2018). These elements included As, Fe, Mg, Mn, and Mo. Additionally, Sn and W were removed because their concentrations were consistently below detection limits. The remaining elements (Sb, Ba, Be, Bi, B, Cd, Ca, Cr, Co, Cu, Pb, Li, Ni, P, K, Se, Ag, Na, Sr, S, Tl, Ti, U, V, Zn, Zr) were geochemically normalized to Al, a lithogenic element, to account for influence of variation in grain size caused by a range of flow conditions (Loring 1991;Kersten and Smedes 2002). This was done to avoid characterization of differences in elemental concentrations that are due to shifting flow energy conditions, independent of source. For example, a single river may supply both coarse-and finegrained sediment of differing elemental concentrations to lakes, and this can confound discrimination of sediment sources. Geochemical normalization is commonly applied to evaluate metals concentrations in lake and river sediment, and Al has been used as the normalizer in this region (Cooke et al. 2017;Kay et al. 2020;Klemt et al. 2020;Owca et al. 2020). We used information in Hugenholtz et al. (2009) to retrieve riverbank sediment samples deposited well before industrial development along the Peace River and Athabasca River. The river-bottom sediment samples and lake surface sediment samples show no indication of enrichment above preindustrial (pre-1920) concentrations (Wiklund et al. 2014;Owca et al. 2020), which allows an approach based, in large part, on contemporary samples to be applied for reconstruction of conditions prior to, and since, upstream industrial development within the river watersheds. LDA was performed using default settings in IBM SPSS Statistics version 27 and forward selection (α = 0.05) based on Mahalanobis distance between the groups to determine a subset of elements that distinguish the three end-members. Leave-one-out crossvalidation was then used to assess the ability of the LDA to categorize sediment samples to the end-members based on elemental concentrations.

Construction, application, and assessment of the BMM
We selected MixSIAR (Stock and Semmens 2016; R Core Team 2020) as the BMM to quantify the source of sediment based on elemental composition, because it provides the flexibility to tailor the model to our specific data set Stock et al. 2018). We used an uninformative prior (1, 1, 1) and the elemental concentrations were entered as the geochemical normalized means and SDs for each of the three end-members (as determined using LDA). Following recommendations by Wynants et al. (2020), we set the Markov Chain Monte Carlo parameters to: chain length = 3,000,000, burn = 2,700,000, thin = 500, and chains = 3. The model was run without fixed or continuous effects. Using the Gelman-Rubin diagnostic and criteria based on Gelman et al. (2013) and Stock et al. (2018), chain length did not require adjustment because < 5% of total variables exceeded 1.05. "Residual" and "process" errors were considered in our model to account for uncertainties introduced by the variable and nonconstant (i.e., pulsed) nature of sediment delivery to lakes across the PAD.
The MixSIAR model was applied to all contiguous freezedried 1-cm thick samples in a sediment core from PAD 30 (length: 45 cm; collected in 2016), and separately again for a core from PAD 31 (length: 38 cm; collected in 2010; Glew 1989), which were analyzed for elemental concentrations at ALS Canada (Waterloo) using the same analytical methods as for end-member samples. The sediment cores were dated using gamma-ray spectrometric measurement of 210 Pb and 226 Ra activities and the Constant Rate of Supply model, as described in Kay et al. (2019). Results were used to assess the model's ability to reconstruct past changes in sediment sources associated with three previously determined, distinct hydrological phases during the past 300 years. PAD 30 (local name: Mamawi Pond; 58.509044 , À111.517280 ) and PAD 31 (local name: Johnny Cabin Pond; 58.495582 , À111.518152 ) are located $ 1 km apart and on opposites sides of Mamawi Creek in the southern Athabasca sector of the PAD (Fig. 1), and thus approximate replicate sites. However, PAD 31 lies at slightly lower elevation than PAD 30, and is slightly more flood-prone (Kay et al. 2019). During the three hydrological phases, mineral matter content in sediment cores varied from $ 40% to $ 90% (Wolfe et al. 2008b;Hall et al. 2012;Kay et al. 2019), which encompasses much of the range in lake sediment composition across the delta. During the first phase, from prior to $ 1700 until $ 1940, the lakes were inundated under a high-stand of Lake Athabasca supported by elevated flow of the Athabasca River. After $ 1940, Lake Athabasca levels declined and both lakes entered a second isolated, closed-drainage hydrological phase of infrequent or absent flooding despite their location adjacent to Mamawi Creek, which became a relict (nonflowing) channel during this period. After 1982, both lakes entered a third hydrological phase when they became increasingly prone to river flooding following an avulsion of the Embarras River into Cree and Mamawi creeks (an event known as the Embarras Breakthrough). This has directed increasing Athabasca River flow via the Embarras River to Cree and Mamawi creeks and to adjacent low-lying areas where these lakes are located. Thus, elemental analysis on sediment cores from PAD 30 and PAD 31 and application of MixSIAR allows us to assess for accuracy and reproducibility of the methodology to reconstruct periods with flooding (under both continuous diffuse flow with occasional pulse flood events pre-1940, and frequent episodic flood events after 1982) and without flooding ($ 1940-1982), and test if the methodology determines that the proximal Athabasca River is the major source of floodsupplied sediment.

Assessment
LDA identified 10 significant elements (Ag, Be, Ca, Li, Pb, S, Sb, Sr, Tl, V) that resulted in clear separation of sediment samples representing the three a priori end-members (Fig. 2). The first axis captured 64.2% of the variation and separated elemental concentrations of sediment supplied by the Peace River, positioned to the right, from sediment supplied by the Athabasca River, positioned to the left. The first axis is most strongly associated with variation in concentrations of Sb, Li, V, and Be. Sb, Li, and V are relatively high in sediment supplied by the Peace River, whereas Be is relatively high in sediment supplied by the Athabasca River. The second axis captured 35.8% of the variation and separated sediment of flooded lakes, positioned low on axis 2, from nonflooded lakes (i.e., local catchment), positioned high on axis 2. Sediment of flooded lakes is relatively high in Pb and Ca, whereas Sr is relatively high in nonflooded lakes that received their sediment from the local catchment. Sample scores of riverbank and river-bottom sediments mostly plotted near the group centroids, as expected given that the Peace River and Athabasca River are the sole source of sediment in these samples. Sample scores of the lakes flooded by the Athabasca River are positioned adjacent to the river sediment samples collected from the Athabasca River, and sample scores of the lakes flooded by the Peace River are positioned adjacent to the river sediment samples collected from the Peace River. Notably, there is no overlap between sediment supplied by the Athabasca River vs. Peace River, and overlap is minimal between surficial sediment of flooded and nonflooded lakes. Overall, use of LDA identifies (1) the source of flood-supplied sediment (Peace River vs. Athabasca River) is the principal gradient of variation in elemental concentrations, (2) variation due to presence or absence of flooding is a secondary gradient, and (3) centroids of the three end-members are distinct.
Cross-validation of the LDA resulted in 90.2% of all sediment samples correctly categorized to their a priori endmember (Table 1). Classification accuracy was slightly higher for the Athabasca River and Peace River end-members (90.0% for Athabasca River, 93.1% for Peace River) than the local catchment end-member (87.9%), consistent with the small amount overlap between sample scores in the LDA from the local catchment end-member with the two river end-members (Fig. 2). Lower classification accuracy of the catchment endmember is not surprising given the lakes are situated in both sectors of the delta. However, the means and sample SDs of the sedimentary elemental concentrations for each endmember are used as input into MixSIAR (Table 2), which are well separated on the LDA as captured by the centroids (Fig. 2). Based on these outcomes, application of the mixing model is likely to provide a robust approach for reconstructing past variation in sediment delivered from the three main sources to lakes in the PAD.
The MixSIAR model was applied to elemental concentrations in stratigraphic records from PAD 30 and PAD 31 spanning the past $ 300 yr to infer temporal patterns of variation in sediment source contribution (Fig. 3). For independent validation, these temporal patterns were compared to variation in organic matter (OM) content, which varies inversely with magnitude of river floodwater influence (Kay et al. 2019). Results from the model produced remarkably similar temporal patterns of variation in sediment proportions at both lakes. During $ 1680-1940, the MixSIAR model indicates that inputs from the Athabasca River dominated sediment composition at PAD 30 and PAD 31 ($ 80%), consistent with relatively low OM content during a period of known prolonged inundation of both lakes under a high-stand of Lake Athabasca. Interestingly, sediment attributed to the Peace River increased to $ 20% during two episodes of this early phase at PAD 30 ($ 1720-1760 and $ 1825-1845). At $ 1940, results from MixSIAR suggest contribution of sediment from the Athabasca River decreased abruptly and supply from the local catchment rose to $ 70-80% at both lakes. These changes coincided with marked rise in sediment OM content and identify onset of the second hydrological phase of PAD 30 and PAD 31 when water levels on Lake Athabasca declined, which caused both lakes to become hydrologically closed due to greatly reduced input of Athabasca River floodwater (Kay et al. 2019). A marked shift in sediment sources is inferred for both lakes in the 1970s, when contribution from the local catchment declined to near-zero and input from the Peace River and Athabasca River increased to $ 57-69% and $ 21-32%, respectively, at PAD 30 and PAD 31. This change was rapid and coincides with the largest ice-jam flood in the instrumental record in 1974, when widespread river flooding inundated the entire delta (Peters et al. 2006). Notably, this flood event is not captured in the OM stratigraphic records. Peace River and Athabasca River sediment continued to be the dominant sources to PAD 30 and PAD 31 until the early 1980s, possibly due to subsequent runoff that transported flood-derived sediment from their local catchments to the coring location. After 1982, sediment from the Athabasca River is inferred to rise suddenly to > 80% in both lakes, coincident with the Embarras Breakthrough which increased frequency of flooding at PAD 30 and PAD 31 by the Athabasca River (Wolfe et al. 2008b;Kay et al. 2019). This shift to the third hydrological phase is also captured by marked decline of sediment OM content due to greater influx of suspended mineral-rich sediment in Athabasca River floodwaters (Kay et al. 2019).

Discussion
Paleolimnological analyses of multiple variables in sediment cores from aquatic basins have provided temporal perspectives about frequency and magnitude of flood events and other hydrological processes across floodplain landscapes (Wolfe et al. 2008b;Lintern et al. 2016). Yet, it remains challenging to integrate results at a landscape-scale because spatiotemporal variations in sediment composition in these dynamic aquatic environments necessitate use of different paleolimnological analyses with varying quantitative reconstruction potential. Thus, development and application of a single analysis that can be universally applied to quantify spatiotemporal patterns of variation in hydrological conditions across floodplains would provide an important advance. Here, we have demonstrated an effective methodology to quantify relative contributions of sediment from two rivers and from local catchments of lakes. We used measurements of elemental concentrations to trace inputs from three different sediment sources and employed a BMM, MixSIAR, to reconstruct variation in hydrological conditions from two lakes with previously described hydrologic history. A distinct advantage of the methodology is that elemental concentrations can be measured across the full spectrum of mineral-rich to organicrich sediment typical of broad gradients of basin hydrology and river connectivity in floodplains.
A fundamental criterion for application of mixing models is characterization of distinct end-members (Smith and Blake 2014;Smith et al. 2015). We relied on geographic location, lake water isotope composition and lake water chemistry data, and knowledge accrued during 20 yr of fieldwork in the PAD, to define the three end-member allochthonous sources of sediment to lakes (Athabasca River, Peace River, local catchment). Clear separation of Al-normalized elemental concentrations among end-members in the LDA resulted in very high classification accuracy (90.2%) by cross-validation due to distinct geochemical fingerprints of Peace River sediment vs. Athabasca River sediment. Such degree of separation of end-member centroids is uncommon among sediment-based mixing models (Upadhayay et al. 2018;Wynants et al. 2020). We achieved this outcome through a substantial effort over many years of fieldwork and knowledge of this complex landscape, yet we recognized this opportunity only after analysis and interpretation of lake surface sediment collected in 2017 for other purposes (Owca et al. 2020). More strategic and efficient sampling is possible elsewhere, which can improve transferability of the approach for use by others with lower investment.
Results of the MixSIAR model applied at PAD 30 and PAD 31 align remarkably well with the three known hydrological phases experienced at these lakes during the past $ 300 yr (Wolfe et al. 2008b;Kay et al. 2019), which capture most of the range of sediment composition for lakes across the PAD. Outputs from MixSIAR inferred that Athabasca River floodwaters were the dominant source of sediment to both PAD 30 and PAD 31 during the past 300 yr, except for a brief period of reduced river influence between $ 1940 and 1982. A novel inference, however, is that the largest magnitude flood on record in 1974 (Peters et al. 2006) resulted in delivery of a substantial proportion of sediment from Peace River floodwaters V 3.08 Â 10 À3 , 3.10 Â 10 À4 3.69 Â 10 À3 , 6.46 Â 10 À4 3.90 Â 10 À3 , 2.95 Â 10 À3 to both lakes. PAD 30 and PAD 31 are proximal (2.7-3.7 km) to the Athabasca River and distal ($ 45 km) to the Peace River, and these results suggest substantial southward penetration of Peace River floodwaters during the 1974 flood which is entirely plausible, based on aerial images of the flood event (Peters et al. 2006). Our methodological framework and application of MixSIAR clearly identifies the most proximal river (i.e., the Athabasca River and its distributary, the Embarras River, following the Embarras Breakthrough in 1982) as the main source of flood-derived sediment, except for the extreme flood of 1974 and possibly also unusually high-magnitude flooding that occurred in the mid-1700s and mid-1800s. Thus, the methodological approach appears to have detected previously unrecognized importance of delivery of floodwater and sediment from the Peace River to the central and southern regions of the PAD during some extreme flood events. Unexpectedly, all three of these intervals of inferred increase of Peace River supplied sediment coincide with rise of sediment OM content, which could occur if the flux of Peace River sediment was insufficient to dilute subsequent within-lake production of organic sediment, and high energy flooding from the proximal Athabasca River did not occur. Overall, for the period of time captured in cores from both PAD 30 and PAD 31 (post ca. mid-1800s), application of the MixSIAR model to elemental concentrations resulted in very similar stratigraphic patterns of inferred sediment proportions from the three end-members at both lakes. Given the comparable hydrological setting of both lakes, we contend the methodology leads to reproducible results that match known temporal shifts in their basin hydrology, yet also identified extended reach of Peace River floodwaters during exceptionally large delta-wide flood events.

Comments and recommendations
For application of the methodology to other floodplains, we recommend that sampling includes river sediment (riverbank, bottom, and suspended sediment) and lake surficial sediment. The river sediment samples are important to anchor the elemental composition of the source material in the mixing model, while the lake surficial sediment samples provide analogs for application to the lake sediment cores. We suggest caution should be exercised when interpreting results from the mixing model to stratigraphic profiles from floodplain lakes because sediment supplied by flood events are deposited on the landscape, which may be remobilized by snowmelt and rainfall events, and lead to inference of river contribution during nonflood intervals. This mechanism could account for inferred substantial Peace River sediment influence in a few samples deposited after the 1974 flood at PAD 30 and PAD 31.
Our methodology advances the important contributions of Wynants et al. (2020) in a new direction. They performed measurements of elemental concentrations in river and lake sediment and applied a BMM to quantify sediment delivery from several tributary watersheds to a single large lake. In contrast, our methodology is capable of quantifying the proportion of sediment contributed to multiple lakes across a floodplain landscape from two rivers during pulse flood events  31 (c, d). The gray-shaded box identifies a closed-drainage hydrological phase of the lakes (ca. 1940-1982) when flooding was rare. A prior hydrological phase is defined by inundation by a highstand of Lake Athabasca and a later hydrological phase occurred as a result of the Embarras Breakthrough when the lakes became prone to frequent flooding from the Athabasca River. OM content and sediment core chronologies for PAD 30 and PAD 31 are reported in Wolfe et al. (2008a) and Kay et al. (2019). and during long-duration inundation, as well as from the local catchment during nonflood intervals. The approach is readily expandable to floodplains with more than two rivers if elemental concentrations of the fluvial sediment differ sufficiently, or can be reduced to those with a single river source. Thus, we anticipate the methodology is transferable to many other floodplains, including the Mississippi River, Tigris and Euphrates rivers, and Yellow River, where flood regimes respond to variation in natural and anthropogenic processes (Richardson et al. 2005;Richardson and Hussain 2006;Nittrouer and Viparelli 2014;Zheng et al. 2017), and where the range of sediment composition may limit application of other paleolimnological measurements.