River Sediment Geochemistry as a Conservative Mixture of Source Regions: Observations and Predictions From the Cairngorms, UK

The elemental composition of sediments in rivers is the product of physical and chemical erosion of rocks, which is then transported across drainage networks. A corollary is that fluvial sedimentary geochemistry can be used to understand geologic, climatic, and geomorphic processes. Here, we predict elemental compositions of river sediments using drainage networks extracted from digital elevation data and erosional models. The Geochemical Baseline Survey of the Environment was used to quantify substrate (i.e., source region) chemistry. Sedimentary compositions in rivers downstream are predicted by formally integrating eroding substrates with respect to distance downstream. Different erosional models, including the Stream Power model and uniform incision rates, are tested. Predictions are tested using a new suite of compositions obtained from fine grained (<150 μm) sediments at 67 sites along the Spey, Dee, Don, Deveron, and Tay rivers, Cairngorms, UK. Results show that sedimentary geochemistry can be predicted using simple models that include the topography of drainage networks and substrate compositions as input. The concentration of numerous elements including Magnesium, Rubidium, Uranium, Potassium, Calcium, Strontium, and Beryllium can be accurately predicted using this simple approach. Predictions are insensitive to the choice of erosional model, which we suggest is a consequence of broadly homogeneous rates of erosion throughout the study area. Principal component analysis of the river geochemical data suggests that the composition of most Cairngorms river sediments can be explained by mafic/felsic provenance and conservative mixing downstream. These results suggest that the elemental composition of river sediments can be accurately predicted using simple erosional models and digital elevation data.

hydrodynamic sorting (e.g., Bouchez et al., 2011Bouchez et al., , 2012Lupker et al., 2016). Consequently, river sediment compositions are used to investigate controls on chemical and physical weathering rates (e.g., climate change, tectonics, and geomorphic processes; Gaillardet et al., 1999;Riebe et al., 2003;von Blanckenburg & Wittmann, 2012). Geochemical studies of fluvial sediment commonly make use of samples along rivers to qualitatively infer climatic or erosional processes upstream. In such studies, assumptions about how the upstream signal is integrated downstream are common, e.g., "let nature do the averaging" (Garzanti et al., 2018;Romans et al., 2016;von Blanckenburg, 2005;Weltje, 2012, and references therein). In this study, we develop quantitative methodologies to objectively test such assumptions using geochemical data.
Predictive modeling of surface processes has become increasingly tractable due to the availability of high resolution topographic data and computationally efficient landscape evolution models (e.g., Braun & Willett, 2013;Hobley et al., 2017;Salles, 2019;Salles et al., 2018). These advancements have made it possible to make testable predictions about a range of observable characteristics such as exhumation rate and fluvial sedimentary flux (e.g., Fernandes et al., 2019). Extending such predictions to geochemical data is a reasonable next step.
The central aim of this study is to develop and test methodologies that can predict the elemental composition of river bedload sediments at any position in a drainage basin. Geostatistical approaches to interpolate geochemical data along stream networks have been proposed (e.g., Kim et al., 2017). Instead we take a deterministic modeling approach, which incorporates assumptions about how physical and chemical signals propagate through fluvial systems. Recently, landscape evolution models have been modified to make theoretical predictions about sediment provenance (e.g., Sharman et al., 2019). Predictions from such models could be tested using, for example, geochemical observations. Such predictions of sediment chemistry along drainage networks are desirable but require information about the chemistry in source regions. In this study, we therefore make use of high resolution regional geochemical surveys, generally produced for environmental monitoring and mineral exploration purposes (Garrett et al., 2008). These inventories can span large regions and contain samples acquired at regular intervals with a consistent methodology, which make them ideal for testing predictive models (Caritat & Cooper, 2016;D. B. Smith et al., 2013).
We begin by describing how the geochemistry of sediment source regions can be constrained using stream-sediment geochemical data. We describe how sediment samples were gathered from higher order rivers to test model predictions. Next we describe a methodology to predict higher order bed material chemistry from topographic data and source region geochemistry. The fitness of this model and its input sensitivity is evaluated by comparing predictions to observations. Finally, we discuss how this approach allows us to simplify the description of sediment sources in our studied region.

Observed Source Region Chemistry: G-BASE Inventory
To predict the composition of bed material in higher order rivers, we need to know the composition of material entering the drainage network. One way to constrain composition of source regions is by interpolating bedrock compositions. A drawback of this approach is that bedrock compositions might not have been affected by chemical weathering. As a result, predictions made using bedrock compositions probably overestimate the abundance of mobile elements. Moreover, sampling bedrock geochemistry at high resolution is challenging. A second way to define source region chemistry is to use the composition of soils. However, because top-and bottom-soil compositions are different, it is unclear which is representative of the material entering the drainage network. The ideal data set for calibrating source region geochemistry would capture geochemical changes in the underlying bedrock as well as incorporate the effect of weathering. A data set which fulfills both of these requirements are geochemical surveys of first-order stream sediments (i.e., streams with very small catchment areas). First-order stream sediments have already experienced weathering, unlike most bedrock, but are not internally heterogeneous like soils. First-order stream sediment geochemistry is also strongly controlled by the geochemistry of the underlying lithology. This LIPP ET AL. 10.1029/2020JF005700 2 of 21 strong relationship between bedrock and stream sediment geochemistry is widely documented (e.g., Caritat & Cooper, 2016;Kirkwood et al., 2016b;D. B. Smith et al., 2013).
We use the high-resolution stream sediment Geochemical Baseline Survey of the Environment (G-BASE) to define substrate (i.e., source region) chemistry. G-BASE was produced by the British Geological Survey to map mineral resources and define environmental baselines. The survey acquired elemental chemistry of sediment samples acquired from low-order streams across the United Kingdom from the 1960s to 2014 (Johnson et al., 2005). The G-BASE survey is distinguished by its high spatial sampling density of one site per ∼2 km 2 . We note that other large scale base-line surveys exist, for example, the National Geochemical Survey of Australia, which samples at a density of one site per 5,500 km 2 (Caritat & Cooper, 2016). High-resolution geochemical surveys such as G-BASE are considered to primarily reflect changes in the geochemical composition of the underlying geology (Everett et al., 2019;Kirkwood et al., 2016b).
The G-BASE survey sampled the fine (<150 μm) fraction of bed material in low-order streams across the UK; sample density is ∼2 km 2 in our study area. A broad suite of elements were measured from each sample. For the remainder of this study, we focus on the following 22 elements, which were consistently recorded in the G-BASE inventory in the studied area: Ba, Be, Ca, Co, Cr, Cu, Fe, K, La, Li, Mg, Mn, Ni, Pb, Rb, Sr, Ti, U, V, Y, Zn, and Zr. For all studied elements apart from Uranium, Direct Reading Optical Emission Spectroscopy was used for the analysis. Delayed Neutron Activation was used to measure Uranium. Full details of the G-BASE sampling, analytical and quality control procedures are given by Johnson et al. (2018aJohnson et al. ( , 2018b. A description of the sampling procedure utilized by the G-BASE survey is given below. The G-BASE sample sites for our study area are shown in Figure 1c. The studied region was sampled by G-BASE between 1977 and1980. It is unlikely that substrate geochemistry has changed dramatically since then for the following reasons. First, the primary control on stream-sediment geochemistry is the underlying lithology which we can regard as constant on the timescales we consider. Second, step changes in geochemical concentrations are not observed between surveys acquired years apart, which suggests continuity of spatial geochemical distributions (Everett et al., 2019). Finally, as far as we are aware, there has been no significant change in the drainage network (e.g., major hydroelectric installations) since the surveys were performed.
The concentration of most elements in the G-BASE data set have a log-normal distribution. Therefore, to generate continuous maps of elements, we interpolated the log 10 transform of the wt% of each element (where wt% = weight percent). The interpolated grids have a resolution of 200 × 200 m (Figures 3e and 3f).
The observations were interpolated using the continuous curvature splines methodology of W. H. F. Smith and Wessel (1990) with tension factor of 0.25. Varying the tension factor between 0 and 1 made negligible difference to predictions.

Study Region
The Cairngorms region was chosen as it contains some of the highest topography in the United Kingdom. It also contains variable lithology (Figure 1b). For example, the Cairngorms massif contains Paleozoic felsic intrusions and mafic igneous rocks intruded into sedimentary and metasedimentary rocks. The diverse array of lithologies, combined with high natural sediment supply provides an opportunity to explore the roles that substrate geochemistry and erosional processes play in determining fluvial sediment geochemistry.
The region was covered by ice sheets during the last glaciation, which receded by 11.5 ka resulting in a layer of till blanketing the region. Reworking of glacial deposits is a significant source of sediment in the region (Ballantyne, 2008). The studied rivers all drain the glacially dissected granite plateau of the Cairngorms massif. The study area has a Warm Temperate/Humid climate (CFb) in the Köppen climate classification (Kottek et al., 2006). A significant part of our studied region is contained within the protected Cairngorms National Park. Outside this region the dominant land use is agricultural. There are no major urban areas upstream of any of the sample sites, although one sample was taken on the river Tay inside the city of Perth and another was acquired from the Don on the periphery of Aberdeen (Localities 55 and 34 in Figure 1d).

Sampling
The G-BASE data set does not contain samples from higher order rivers and hence cannot be used to evaluate our model or investigate the length-scales at which downstream bed load geochemistry varies. Consequently, a new inventory of 67 samples was acquired for this study along higher order rivers draining eastern Scotland, UK, in August 2019. We focus on the Spey, Tay, Dee, Don, and Deveron rivers, which drain the Cairngorms massif and surrounding region of northeastern Scotland, UK ( Figure 1).
Sample sites in the Spey, Tay, Dee, Don, and Deveron catchments were chosen to include the main trunk from the upper reaches to the mouth and most major tributaries. Where possible, samples were gathered before and after a major confluence from both channels. The exact location of sample sites was determined by the practical need to access the river. To maintain consistency between the model input data and our test LIPP ET AL.  data set, the sampling procedure was identical to the original G-BASE approach and used the same equipment ( Figure 2a). The standard sampling procedure used for the G-BASE study, and repeated by this study, is described below.
We sampled the 150 Pm fraction of bed-material. Bed-material was extracted by shovel from the active channel of the river. Examples of sample sites are shown in Figures 2b-2d. For practical reasons, samples were gathered from bed-material close to the bank, within ∼2 m (e.g., Figures 2c and 2d). If multiple shovel loads were required, as much as was possible, they were extracted from the same point in the channel bottom. The extracted sediment was first shaken through a 2 cm metal grill to remove pebbles. Subsequently, the wet sediment was deposited on a sieve stack on top of a fiberglass collecting pan. Using rubber gloves, the sediment was rubbed through a 2 mm nylon mesh and washed through with a small amount of river water. This upper mesh was removed and the <2 mm sediment fraction was then rubbed and washed through a 150 μm nylon mesh into the collecting pan beneath. This procedure was repeated until ∼100 g of sediment had collected in the pan. The collecting pan was left undisturbed and covered for 20 min to allow sediment to fall out of suspension. Excess water, which contained some very fine sediment in LIPP ET AL.  suspension, was carefully decanted leaving a slurry of sediment behind. This slurry was homogenized and poured through a funnel into a labeled reinforced paper sample bag. Any sediment residue was then washed into the sample bag with a small amount of water. The bag was sealed and placed in a sealed plastic bag. A video showing the sampling procedure can be found in the Supporting Information.
To limit cross-contamination, the sampling kit was washed in river water, downstream of the site, before and after sampling. Sample numbers were pre-allocated in a randomized order. For example, localities 1, 2, and 3 had pre-allocated sample numbers CG020, CG062, and CG044, respectively. This randomization was performed to avoid systematic bias during preparation and analysis, when samples were handled in numeric order (Johnson et al., 2018a). At the end of each field day, samples were dried in the paper bags to remove excess water until they had the consistency of modeling clay. Special care was taken to avoid LIPP ET AL.
10.1029/2020JF005700 6 of 21 contaminating samples at this stage. The samples, still in the paper bags, were returned to their polythene bags and resealed. After the sample was concluded, samples were freeze-dried to remove any remaining moisture before short-term storage prior to analysis.

Geochemical Analysis
The freeze dried samples were disaggregated with a rubber mallet and homogenized by cone-quartering. Twenty grams of the homogenized sample were then powderized in an agate ball mill. Approximately 0.25 g of the milled powder was extracted and then accurately weighed in Savillex tubes for digestion. The samples were predigested in HNO 3 to remove organic compounds. Subsequently the samples underwent hotplate mixed acid (HF, HNO 3 , and HClO 4 ) digestion over a ramped heating procedure. After digestion the samples were resolublized in aqueous HNO 3 and H 2 O 2 . The liquid samples were then analyzed for a full suite of elements by Inductively Coupled Plasma Mass Spectrometry. Comparison to external and internal standards indicates that, as expected for our method, some elements (Zr, Y, and Ti) reported slightly lower concentrations in our standards than expected suggesting a small residue of particularly chemically resistant minerals such as Zircon. The majority of other elements in the standards were successfully reproduced. All sample preparation and analysis was carried out at the BGS Inorganic Geochemistry Facility in Keyworth, UK. The new data set of chemistry for higher order rivers is provided (see Acknowledgments for access details).

Nested Duplicate and Replicate Analysis
The goal of our study is to identify and predict geochemical differences between sample sites, that is, intersite variability. These intersite differences may however be obfuscated if there is geochemical heterogeneity within each sample site, or introduced by our analytic procedure. To investigate the ability of our procedure to identify regional geochemical signals, we conducted a nested sampling procedure ( Figure 2e). Nested sampling procedures are commonly used in geochemical surveying to quantity the amount of intersite, intrasite, and intrasample variability and hence evaluate the success of the sampling and analytical procedure.
We implemented this nested procedure by gathering duplicate samples at four randomly chosen localities ( Figure 1d). The duplicate samples were gathered exactly as described above but at a distance of ∼100 m from the previous sampling point. Prior to powderization, each homogenized duplicate sample was then split into two "replicate" samples to investigate intrasample variability such that each duplicated site yields four analyses (Figure 2e).
We use a nested analysis of variance (ANOVA) to quantify the intersite, intrasite, and intrasample variance for each of the studied elements in the 16 replicate/duplicate samples (Garrett, 1969). Nested ANOVA partitions an element's variance (the spread of values for a given element, given by  2 ) into specified hierarchies (e.g., Figure 2e). In doing so, it reveals the contribution from intersite, intrasite, and intrasample variability toward an element's variance. Because the elements had a log-normal distribution, we perform ANOVA on each element after applying a log 10 transform.

Elemental Compositions
The downstream composition of river sediment can be predicted using information about the properties of eroding substrate (e.g., elemental composition of substrate) and incision rates,   / z t. In this study, the former is quantified by the G-BASE data set. A suite of observational and theoretical approaches exist to calculate incision rates and sedimentary fluxes at a range of spatial and temporal scales (e.g., Holbrook & Wanas, 2014;Stephenson et al., 2014;Syvitski & Milliman, 2007). For example, cosmogenic dating of fluvial terraces and erosion of radiometrically dated basalt flows constrain incision rates at spot locations on time scales of 1 kyr-1 Myr (e.g., Karlstrom et al., 2008;Stucky de Quay et al., 2019). Erosional models are also frequently used to calculate incision rates continuously across a landscape.
First, for simplicity, we consider some spatially variable and measurable geochemical characteristic of sediment, C (e.g., wt% MgO) in a single river. By assuming instantaneous transport of eroded sediment, the composition of river sediments C sed along a river of length L is given by This expression calculates the contribution of source regions to the value of C sed -by integrating the product of composition and incision rate with respect to distance downstream between positions x and L. The contribution of elements from each eroding patch of the landscape are then normalized by total sedimentary flux.
It is straightforward to generalize this scheme to two (spatial) dimensions by integrating all of the upstream flowpaths from the point x. In this way, the composition of river sediments can be predicted using all of the eroding catchment upstream. To calculate a cumulative volumetric sedimentary flux, Q sed , we integrate the incision rate across the entire upstream area, A, from a point, that is, (2) We make use of topographic data, z(x, y), and continuous maps of source region geochemistry, C(x, y), to parameterize these equations. This approach integrates all upstream regions that could contribute to compositions downstream. The two-dimensional extension of Equation 1 is therefore: For this study, we choose to use erosional models to calculate incision rates due to our requirement for incision to be defined at all points in a landscape. Sparsely located point observations would be inappropriate for this purpose.

Erosional Models
We first examine predictions from the widely used Stream Power erosional model, which relates incision rates to discharge and slope,   / z x, through a power law relationship (e.g., Howard & Kerby, 1983;. This model is simple, commonly used and can be readily extended to consider, for example, "diffusive" geomorphic processes. Such diffusive processes can modify hillslope shapes in addition to the dominant fluvial "advective" processes (e.g., Rosenbloom & Anderson, 1994). Using upstream drainage area, A(x), as a proxy for river discharge results in the following formulation of fluvial incision rate along a single river where v, m, and n are parameters that can be calibrated using the topology of drainage networks and independent geologic data (e.g., Quye-Sawyer et al., 2020;Stock & Montgomery, 1999). This formulation is often presented as E(x,t) = −vA m S n , where S is the slope and E is the erosion rate. The parameter v is often interpreted as representing bedrock erodibility, which could vary as a function of lithology, structural weaknesses, and climate (Howard & Kerby, 1983;Whipple & Tucker, 1999). The term A(x) m is a proxy for the power exerted on the river bed by the channel, where the value of m is controlled by basin hydrology. Specifically, m depends on both how channel width scales with discharge and how discharge scales with upstream area. The parameter n controls how the knickzone retreat rate depends on slope. For example, n > 1 implies that steeper knickzones migrate faster than shallower ones for the same upstream area (Pritchard et al., 2009).
To calculate erosion rates across an entire two-dimensional landscape we solve a more general form of the Stream Power model, such that, where z is the steepest cell-to-cell descent (see e.g., Salles et al., 2018).
LIPP ET AL.

10.1029/2020JF005700
Inverse modeling of large inventories of river profiles has shown that n ≃ 1 and 0.2 ≤ m ≤ 0.7 produce theoretical river profiles that best fit observed longitudinal river profiles in Eurasia, North and South America, Africa, and Australia (e.g., Fernandes et al., 2019;Paul et al., 2014;Rudge et al., 2015). Therefore, we start by assuming spatially constant values of n = 1 and m = 0.35 (Figures 3a-3d). We start by assuming that v is also constant; as such the predicted composition of sediment is independent of its specific value. However, as an intermediate step in our calculations v must be specified so we, arbitrarily, use a value of 3.62 m 1−2m Myr −1 following Paul et al. (2014). We acknowledge that a range of different values of m and n are preferred (see Royden & Perron, 2013;, and references therein). As a result, we explore the sensitivity of our model by systematically varying erosional parameter values (e.g., n = 2, m = 0.35; n = 1, m = 0.6). We also consider a simple scenario in which incision rates are constant across the landscape. The appropriateness of chosen parameters is assessed by calculating misfit between observed and predicted chemical compositions along higher order rivers.

Implementation
Equations 3 and 5 are implemented using the LandLab package in python 3 (Hobley et al., 2017;Van Rossum & Darke, 2009). The SRTM1s topographic data set was gridded to 200 × 200 m squares following a cylindrical equal-area projection centered on our study area using GMT 5.4.5 (Farr et al., 2007;Wessel et al., 2013). Performing the calculations with a finer resolution grid significantly increases computation time. Depressions in the topography were filled using the "priority-flood" algorithm and flow directions are computed using the "D8" algorithm (see Barnes et al., 2014). Code to reproduce our implementations of Equations 3 and 5 are provided (see Acknowledgments for access details).

Nested ANOVA
The results of the nested ANOVA showing the proportion of variance in each hierarchy are shown in Figure 2f. The intersite variance accounts for >90% the total variability for most elements. This result indicates that regional geochemical signals dominate over local variability. In most instances, the intrasample variance, which mostly reflects measurement uncertainty in our analytical procedure, is negligible. For Be, the relatively high-intrasample variance is also reflected in the digestion repeats where Be shows variability between the repeats. Given that the intra-sample variability is so small for all elements, we use the arithmetic mean of the two replicates to represent the duplicate value for subsequent analyses. We conclude that our sampling and analytical methodology is successful in identifying regional geochemical signals between sample sites. Figure 4a shows observed sediment Mg concentrations and model predictions. Figure 5 shows results for elements that are well matched by predictions. Figures 6 and 7 show results, alphabetically, for other elements predicted by our model. Figure 8 shows results for elements poorly matched by the model.

Evaluating Predictions
The model predictions show variability in sediment geochemistry, even across relatively small areas. For example, the Don has a higher predicted Mg concentration at its mouth than the Dee, despite being separated by only ∼4 km (Figure 4). River chemistry is predicted to evolve along river channels. The Dee, for example, has very low predicted concentrations of Mg in its upper reaches but sediment becomes progressively more enriched in Mg toward its mouth. Variability on this scale can be observed in the predictions for all elements. This variability is validated by the nested ANOVA results which indicate that the primary geochemical signal in bed material exists between different rivers and topologically distinct points within a drainage network. Visual inspection of maps showing observed and theoretical sediment chemistry indicates that our simple model, which incorporates substrate chemistry and drainage networks, can reliably predict river sediment chemistry for many elements.
We formally compare the observed and predicted elemental concentrations by calculating root mean square (RMS) misfit and R 2 values. These statistics are calculated using the log 10 transformed data for two reasons. First, the data have a log-normal distribution, and second, if raw data (wt%) were used, strong heteroscedasticity was observed, which violates the underlying assumptions of the chosen statistics. If the log-transform was not used, undue weight was applied by the regression algorithm to samples with high concentrations skewing the results. To extract model predictions for each observation, we select the closest model cell that exists in a channel (defined as a cell with upstream area >25 km 2 ). This procedure was LIPP ET AL.
10.1029/2020JF005700 10 of 21 checked manually for each locality. These statistics are given for every element alongside the results maps (Figures 4-8).
The misfit is calculated as the difference between the (log transformed) predicted and observed values, that is, log 10 C o − log 10 C c . In most instances, the misfits of the elements are distributed around 0 indicating limited bias. Ca is the only element which is systematically under-predicted by our model, despite having a high R 2 value (Figure 5d). Eight of the studied elements (K, Mn, U, Zr, Li, Co, Ti, and Y) are overestimated by our model. This over-prediction can be partially explained by our analytical procedure which can leave strongly resistate minerals only partially digested. Zr, Y and Ti are often hosted in resistate minerals such as Zirconium and Titanite and hence could be more abundant in the higher order river samples than our results indicate. Even in these instances of over-or under-prediction the fitted regression lines are parallel LIPP ET AL.
10.1029/2020JF005700 11 of 21 or sub parallel to the target 1:1 line. Exceptions to this are Pb, Mn, and Y. Plotting the misfits spatially shows no spatial bias, such as shown for Mg in Figure 4b. Maps of misfit for all elements can be found in the Supporting Information.
The RMS and R 2 values for every element are displayed graphically in Figure 9a. For multiple elements such as V, Rb, U, and K (in total 9 out of 22, or 41%), the predictions account for >70% of the observed variability ( Figure 5). For the majority of elements, (in total 15 out of 22, or 68%), the R 2 value is > 0.5, which indicates that the model successfully captures the majority of the observed variability.
LIPP ET AL.

Model Sensitivity
We quantified the extent to which our model predictions are sensitive to different erosional parameters by changing the values in the model and observing resultant model predictions. We find that our model is insensitive to different erosional parameter values. Changing the value of the exponent n (but keeping all other parameters the same) between 0.05 and 2 in Equation 4 has a very limited effect on the predicted geochemistry and only very weakly affects the model fit where the maximum change in R 2 , for Ca, is from 0.73 to 0.70 (Figures 4c, 9b, and 9c). Similarly, using a different m value of 0.6 also had a small effect on model fit, where the maximum change of R 2 was for Sr from 0.70 to 0.69 (Figures 9b and 9c). Assuming a completely homogeneous incision rate across the studied area resulted in a small increase in R 2 value for LIPP ET AL.
10.1029/2020JF005700 13 of 21 17 out of 22 of our studied elements (Figures 9b and 9c). The maximum increase in R 2 was observed for Ba from 0.53 to 0.57.
By contrast, our model is very sensitive to the geochemical input. We explored this sensitivity first by spatially randomizing the original G-BASE input prior to interpolation and secondly by using an equally sized rectangle of G-BASE data from an arbitrarily chosen part of the United Kingdom (Wales). These alternative inputs are shown, for the case of Mg, in Figure S1 of the Supporting Information. Using these different geochemical inputs results in much lower R 2 values for all elements with all R 2 values dropping to below 0.2 and an appreciably higher RMS misfit in all instances. Consequently, the success of our model when using the "true" G-BASE data set is not simply coincidental and reflects a successful integration of the upstream geochemistry. This result emphasizes that input geochemistry must be known in some detail to predict downstream geochemistry.

Explaining Successful Prediction of Sediment Geochemistry
The results of the ANOVA indicate that regional (kilometer to tens of kilometer) variations in river sediment geochemistry dominate over local (meters to hundreds of meters) heterogeneity ( Figure 2). While processes such as hydrodynamic sorting and grainsize are generally agreed to affect sediment geochemistry locally, these effects appear to be subordinate to larger regional signals.
The successful model predictions indicate that river sediment geochemistry is primarily controlled by conservative mixing of heterogenous source regions. The lack of bias suggests that in-transit modification of LIPP ET AL.
10.1029/2020JF005700 14 of 21 sediments by physical and chemical weathering through abrasion and attrition is likely a secondary process, consistent with a recent study of Boron isotopes (Ercolani et al., 2019). The elements which are over/ under-predicted by the model, as discussed in the Section 3, could be exceptions to the general conservative behavior. For example, the concentration of Ca is consistently under-predicted by the model. It is possible that cation exchange with the dissolved load could possibly account for this effect, as dissolved Ca is observed to adsorb onto riverine sediments (Cerling et al., 1989). We note that the study area has a relatively temperate climate and does not contain protracted periods of sediment storage in floodplains. Both these factors reduce the possibility of in-transit chemical and physical weathering. The composition of sediments in the Amazon and Ganges rivers, which do have large tropical floodplains, has been interpreted in terms of in-transit sediment modification which could indicate nonconservative behavior (Bouchez et al., 2012;Lupker et al., 2012). Nonetheless, in the absence of protracted sediment storage, we suggest that downstream changes in sediment geochemistry should be interpreted primarily as a result of mixing of different source regions and not in terms of changing intensity of a particular climatic or geomorphic process. While we consider only elemental geochemistry in this study, similar consideration could perhaps also be applied to sediment hosted isotopic proxies (e.g., Lithium isotopes,  7 Li).
LIPP ET AL.

Model Limitations: Grain Size, Temporal Evolution, and Erosional Regimes
In this study, we account for the effect of sorting of sediments by simply extracting the fine grained portion of bed material. The sampling procedure was designed to follow that used to generate the G-BASE data set. However, a drawback of this approach is that grainsize induced artifacts could have been introduced into the compositional data set. Our approach implicitly assumes that the grainsize distribution of the sampled bed material fraction ( 150 Pm) is the same at each sample site. As a result, variations in the distribution of grainsizes across the studied region could impart compositional variability that is not captured by the model. Given that grainsize distributions are generally observed to vary in most riverine bed material, some of the geochemical variability not explained by our model is likely due to this effect. Indeed, some of this variability will contribute to the intrasite variance identified by the ANOVA (Figure 2f). More sophisticated sediment sampling procedures based on hydrodynamics would probably better resolve the larger regional signals from the grainsize induced "noise" (e.g., Baronas et al., 2020;Lupker et al., 2011).
The nested ANOVA was used to quantify the possible role of local (tens to hundreds of meters) geochemical heterogeneity. However, we are not able, in this study, to quantify the possible magnitude of temporal changes in sediment geochemistry at each sample site. For example, localized geomorphic events such as landslides can induce temporary changes in sediment supply from particular parts of a catchment. In addition, seasonal changes in river hydrodynamics (e.g., flood events) may impart a seasonal grainsize distribution variability that would translate into geochemical changes as discussed above. A further concern is that the recent actions of humans have affected sediment dynamics in our study region. Globally, humans have accelerated soil erosion rates via deforestation and the expansion of agriculture and locally, dams and other human structures reduce sediment supply (Syvitski & Milliman, 2007). The studied region contains both large areas of agricultural land-use and also two major dams which could introduce sediment supply distortions not accounted for by our model. The protected status of the Cairngorms National Park could reduce the magnitude of these potential anthropic effects. These temporal and anthropic fluctuations will contribute to the geochemical variability not explained by our model. We note that cosmogenic nuclide studies suggest that the average residence time of sediment in drainage basins is likely to be of the order of hundreds to thousands of years (Repasch et al., 2020). These observations suggest that the composition of sediment in drainage basins might integrate processes operating on timescales of hundreds to thousands of years. Hence, while temporal and anthropic induced geochemical variability likely exists at some scale, it is probably subordinate to the dominant spatial trends predicted by the model and confirmed by measured compositions.
Predicted compositions are very sensitive to geochemical input, which suggests that the first-order control on river sediment geochemistry is the distribution of source rocks in their catchments. The dependence on source composition is also a consequence of drainage basins containing fluvial networks that have specific topologies (i.e., geometries and spatial relationships). In contrast, we found that model predictions were insensitive to the different fluvial erosional models we tested. The linear (n = 1) and nonlinear (n = 2) dependencies of incision on slope fitted the data equally well. Additionally, by changing the value of m we explored sensitivity of predictions to changing basin hydrology scaling relationships. The minor differences in R 2 values suggest a weak dependency on hydrology. The most successful (highest R 2 values) erosional model we apply is homogenous erosion, that is, a constant lowering rate across the landscape. While in this location constant lowering rate generates accurate predictions, it is almost certainly not a reasonable proposition in many settings where erosion is highly variable (e.g., in active orogens). We suggest that these results are a consequence of the choice of study area, which probably had a broadly stable Holocene erosional regime and is not tectonically active. In other regions, where the erosion rates are more spatially variable (e.g., across faulted regions which generate differential uplift or with strong climatic gradients), we might expect to observe divergence between geochemical predictions from different erosional models. An analogous study to our own, but in a more tectonically active area, would likely allow different erosional models to be better discriminated. We also acknowledge that more complicated erosional models could be used, such as incorporating diffusive hillslope processes. However, given that the predicted geochemistry is relatively insensitive to the particular erosional models tested, it is unclear how additional complexity would provide further insight.

Simplifying the Controls on Cairngorm River Sediment Composition
A key question is what controls the spatial distribution of elemental concentrations recorded in river samples. So far we have considered each element separately, which effectively results in a 22-dimensional data set. This high dimensionality makes interpretation challenging. To simplify interpretation, we apply principal component analysis (PCA) to the higher order river data set following a centered log-ratio transformation (Aitchison, 1983). PCA is a widely used dimension reduction technique and is increasingly applied in sedimentary geochemical analyses (e.g., Vermeesch & Garzanti, 2015). PCA works by projecting the raw  data onto a smaller number of principal components along which the variance is maximized. Hence, by selecting components which contain large amounts of variance, the dimensionality of the data set is reduced.
For the higher order river sampling data set, 67% of the total variance is contained on the first principal component (Figure 10d). Inspection of the loadings, that is, the weightings of the original variables, on the first principal component suggests that it corresponds to a discrimination between rocks that might broadly defined as "mafic" and those that might be called "felsic" (Figure 10e). Generally, elements with positive loadings are enriched in mafic rocks (e.g., Cr, Mg, and Ni) whereas elements with negative loadings are enriched in felsic rocks (e.g., Be, K, and U). This distribution is confirmed by projecting the G-BASE data onto this component. The resulting map, shown in Figure 10a, clearly defines different geological units (Figure 1b). The felsic intrusions are highlighted in blue, relative to the more mafic rocks shown in red. We emphasize for this discussion we use "felsic" and "mafic" as two endmembers which exist on a geochemical continuum, not specific lithologies. Hence, the geochemical composition of sedimentary rocks (which span much of our study region) can also be projected based on their geochemistry onto this score.
This result suggests that the Cairngorms river sediments can be well described simply in terms of having either broadly mafic or broadly felsic source regions. Projecting the predicted and observed sediment compositions onto this principal component confirms this result (Figure 10b). The Spey and Dee rivers have negative PCA scores indicating a generally felsic catchment. By contrast, the Don and Deveron have positive PCA scores at the mouth, indicating a generally mafic catchment. This analysis emphasizes the sensitivity of river sediment compositions to the source region geochemistry. The Don and Dee have contrasting mafic and felsic sediments but for much of their length they share a common watershed and are separated at the mouth by only ∼4 km.

Applications and Further Work
Geochemical surveys such as G-BASE are widely used as baselines for environmental monitoring (e.g., Young et al., 2016, and references therein). Higher order rivers are rarely included in sampling campaigns (for an exception, see Fordyce et al., 2004). Consequently, a generic baseline for environmental monitoring of higher order river sediments is lacking. Extending baseline surveys to include sediments contained in higher order river sediments, which, as this study shows, integrate geochemistry of their upstream region, would therefore fill a present gap in baseline applicability. Predictive models such as the one we propose, which utilize survey data as input, could be used to make baseline predictions for higher order rivers.
We have demonstrated that sediment geochemistry can be predicted in rivers from a model that assumes conservative mixing given well constrained source region geochemistry. Here, incision is described by a Stream Power model, although all conservative mixing models tested in this study are similarly successful. This predictable behavior of river sediments suggests that their compositions could be formally inverted to obtain maps of source-region geochemistry in drainage basins.
A further prospect implied by our results is that landscape evolution models could be parameterized using sedimentary geochemical data. In this study, we have only considered a relatively small number of different parameters, changing the values of n and m individually. However, performing a systematic sweep of the n and m parameter space, to identify the parameter set that best fits the higher order sediment geochemistry would be a novel way to parameterize a landscape evolution model. A potentially limiting scenario is the need for high-density geochemical data in the source regions. Nonetheless, many national geochemical surveys exist (e.g., in the US; mrdata.usgs.gov/geochem) and geochemical surveying is frequently performed in areas of potential mineral resources by the mining industry. There is therefore likely a large amount of legacy geochemical data that could be utilized, in combination with specifically gathered higher order river data, to parameterize erosional models.
More generally, our approach can be tested and applied anywhere with suitable upstream geochemical data. While G-BASE is, as far as we aware, the only national geochemical survey of such a high density, the approach we propose could plausibly be adapted to surveys with lower density sampling. The specific sampling density required will, however, be dependent on the expected geochemical variability of the studied region. Where geochemical data are not available, it could be possible to parameterize geo-logical maps, which are more widely available, in terms of geochemical compositions. Finally, recent advances have been made in using machine learning to predict substrate chemistry from a range of more widely available remote sensing and potential field data (e.g., Kirkwood et al., 2016a;Wilford et al., 2016).

Conclusions
We develop a model to predict the geochemical composition of higher order river sediments using erosional models and maps of source region geochemistry. This scheme is tested in a case study of the Cairngorms, UK. Statistical analysis of point measurements gathered across the region show that regional geochemical variability dominates over local heterogeneity. The model is successful in explaining the majority of the observed variance for 15 of our 22 studied elements. For our chosen region, the predicted sediment geochemistry is insensitive to the erosional parameters chosen, which is likely due to the homogeneous climatic and tectonic regime of our study area. By contrast, we find that sediment geochemistry is highly sensitive to the spatial distribution of the source region geochemistry. In the Cairngorms river, sediment composition is primarily set by conservative mixing of heterogeneous source regions, which can be predicted using simple erosional models. In our study region, river sediment geochemistry can be additionally described simply in terms of mixtures of source region geochemistry which exist on a continuum between mafic and felsic endmembers. Our predictive scheme indicates that quantitative modeling, and inversion, of fluvial geochemical signals is a promising avenue of further research.