Q uantifying excess heavy metal concentrations in drainage basins using conservative mixing models

identify excesses


Introduction
River sediment compositions are an important record of natural and anthropogenic processes upstream (e.g., Caracciolo 2020). For example, increased heavy metal concentrations can indicate river reaches polluted by mining (e.g., Sengupta 2021). They can also provide cost-effective means to map the composition of sources areas and are used as tools for exploration and to identify reaches requiring remediation (e.g., Lechler et al. 1997;Luıs et al. 2011;Wijaya et al. 2013). Methods to identify theoretical or end-member sources of material are an active line of inquiry (e.g., Cicchella et al. 2022). The goal of this study is to use river sediment compositions to identify the actual locations of pollutants and other sources of material upstream. The composition of river sediments depends primarily on the mixture of elements in upstream source regions (e.g., Ercolani et al. 2019;Caracciolo 2020;Lipp et al. 2020). Secondary processes such as weathering, cation exchange and sorting play moderating roles, alongside any compositional changes due to anthropogenic input (e.g., Bouchez et al. 2011;Lupker et al. 2012;Tipper et al. 2021). Determining the provenance (i.e., concentrations as a function of location) of river sediment is, however, complicated by mixing of heterogeneous sources downstream. In this pa-per we demonstrate how new computational strategies can be used to 'unmix' sediment compositions to identify locations of natural and anthropogenic sediment sources. We identify three challenges associated with using spot measurements of river sediment chemistry to identify and quantify natural and anthropogenic sediment sources.
First, source regions are chemically heterogeneous and can include natural (e.g., rocks and in situ derivatives) or man-made (e.g., industrial waste) material. Source region heterogeneity means that the composition of sediment in rivers downstream is also likely to be heterogeneous. Consequently, actual source region compositions should be incorporated into mixing models (cf. Luıs et al. 2011;Goswami et al. 2021). Secondly, spot measurements of concentration need to be interpolated to generate continuous predictions across drainage networks. Interpolating observations across drainage networks is not trivial as most methods consider only geographic/cartesian distance, which is inappropriate for drainage networks (see e.g., Kim et al. 2017). Finally, to quantify anthropogenic contributions to the concentrations of elements along rivers we first need to isolate contributions from naturally occurring processes (e.g., mixing of weathering products). Preprint -Quantifying excess heavy metal concentrations in drainage basins using conservative mixing models 2 In this paper we first demonstrate how simple conservative mixing models can be used to generate 'natural' baseline concentrations of elements in rivers using high-resolution geochemical maps of source regions. Second, we show how these baselines can be used to identify excess heavy metal (Pb, Cu, Zn) concentrations downstream. We then use an inverse methodology to explore if such excess concentrations are likely to be of anthropogenic origin. The Upper Clyde Basin, near Glasgow, Scotland is used as our case study (Figure 1).

Strategy
The methodology presented in this paper is summarised in Figure 2 and has three parts. First, we generate predictions of element concentrations along rivers using conservative mixing models that assume the composition of sediment is set by natural 'geologic' input ( Figure 2: Forward Model). Formally, 1185 measured concentrations of first-order stream sediments are integrated with respect to upstream area to predict compositions of river sediments downstream. We assume that these predicted downstream concentrations are 'natural' baselines, with which measured concentrations along rivers can be compared. For some elements, e.g., Mg, Sr, K, observed concentrations are tightly distributed around baseline predictions, thus they are considered to not have significant contributions from anthropogenic sources. Other elements, including the heavy metals Pb, Cu, and Zn, display discrepancies between observed and baseline concentrations. These elements are interpreted to be influenced by anthropogenic sources. The second part of our study is focused on quantifying the excess material in source areas required to match observed concentrations downstream. It involves generating upstream geochemical maps that best-fit the relatively small inventory of 60 concentration measurements along the Clyde river. The best-fitting model is identified by seeking smooth synthetic upstream (source area) concentrations, which are mixed downstream following the forward modelling procedure described above (Figure 2: Inverse Model). Finally, differences between calculated baselines and the concentrations that minimise misfit to observations along the river are used to quantify excess, likely anthropogenic, input. The results are used to identify locations along the river where heavy metal concentrations exceed toxic limits and locations of potential contaminants in sources areas. In the following section we discuss the choice of study area, data used to parameterise the computational models and the mathematics that underpin the methodology.

Study Region
We focus on the upper reaches of the Clyde river, Scotland, where the composition of sediments at 60 localities along the main channel upstream of Glasgow were measured by the British Geological Survey (BGS) as part of the Clyde Urban Super Project (CUSP; Figure 1; Fordyce et al. 2017). The upper Clyde Basin covers an area of 1873 km 2 and is mainly composed of sedimentary rocks: sandstone and wackes (Figure 1b). Felsic and mafic intrusions can be found throughout the region including ∼ 10 km 2 of extrusive igneous rock southwest of Glasgow.
The Upper Clyde region has been home to intensive mining and related industry for the last several centuries (Skillen 1987). Economic decline after World War I saw much of the industry abandoning the city, leaving behind high rates of chronic disease and one of the lowest life expectancies in Europe (Kintrea and Madgin 2019). Environmental remediation was a key aspect of the regions successful revitalisation (Maver 2019). The CUSP, a large scale geochemical surveying campaign to understand the current state of the environment in the Clyde Basin and estuary 2 , was created to support these efforts (Campbell et al. 2010;Appleton et al. 2013;Morrison et al. 2014). In ten years, this campaign gathered 5943 stream, soil and water samples and analysed them for a range of parameters including elemental chemistry .
In this study, river sediment samples collected along the main channel of the Clyde during the CUSP are combined with firstorder stream samples from the BGS's Geochemical Baseline Survey of the Environment (G-BASE) to quantify excess heavy metal concentrations in the region through conservative mixing models ( Figure 2).

G-BASE
Two pre-existing geochemical data sets are used in this study. First, elemental concentrations in stream sediments, recorded in G-BASE, are used to quantify source region composition. G-BASE was a multi-decade  sampling campaign with the goal of generating a geochemical baseline for the UK (Johnson et al. 2005). Sediment samples from first-order streams were obtained with an average sampling density of 2 km 2 . Full sampling and analytical protocols for G-BASE are given in Johnson et al. (2018b) and Johnson et al. (2018a). Figure 1c shows the location of the 1185 samples within the Clyde basin; note that the Glasgow urban area has the lowest sampling density. Through sieving, the grain-size fraction < 150 µm was extracted and analysed using X-ray Fluorescence (XRF) and direct reading optical emission spectroscopy. Figure 3 shows that the elemental concentrations of first-order streams are highly correlated with underlying lithologic units (e.g., high Mg in the south-east coincides with layers of wacke). We consider stream sediment compositions as source region concentrations in this study because chemical weathering has already occurred prior to the sediment reaching the streams.

Upper Clyde Sediments
The second dataset is a suite of sediment elemental concentrations measured at 60 localities along the main Clyde channel, and select tributaries. These samples were collected as part of the Upper River Clyde Sediment Survey, a section of the CUSP. These river bed samples were collected along the main trunk of the River Clyde and at the mouths of major tributaries (Smedley et al. 2017). They have been analysed in accordance with G-BASE protocols (grain-size fraction < 150 µm). Their distribution is shown in Figure 1d. Note that additional CUSP samples were also collected further downstream along the Clyde (Jones et al. 2017). These samples are not considered further Preprint -Quantifying excess heavy metal concentrations in drainage basins using conservative mixing models 4 Preprint -Quantifying excess heavy metal concentrations in drainage basins using conservative mixing models 5 here because a different grainsize fraction (< 2 mm) relative to the rest of the dataset was used.

Numerical Models
Two computational methods are used in this study. First, continuous concentrations of elements along rivers ('natural' baselines) are predicted using a simple conservative mixing (forward) model (Lipp et al. 2020). Secondly, the CUSP (higher-order river) samples are inverted ('unmixed') to predict concentrations of elements in upstream source regions. For both the forward and inverse model a stream network must be defined. First, the SRTM1s DEM is downsampled onto a 100 × 100 m grid and loaded into the Python package LandLab 2.3.0 (Hobley et al. 2017;Barnhart et al. 2020). Secondly, the topography is filled using the priority flood method to avoid artificial depressions (Barnes et al. 2014). Finally, the D8 flow-routing algorithm is used to generate the stream network from the SRTM 1s topographic data (O'Callaghan and Mark 1984;Farr and Kobrick 2000). Rivers in this study are defined by drainage areas > 8 km 2 .

Forward Model
Assuming conservative mixing, there are two controls on the downstream concentration D(x) of an element, upstream source concentration, C(x), and the vertical erosion rate ∂z/∂t, which yields the following relation where x is a location along a streamline, A is the upstream area, z is elevation and t is time. In other words, river sediment compositions are a function of upstream source region compositions mixed downstream. Forward modelling in the Cairngorms, Scotland, suggests that different erosional models have a minor impact on the conservative mixing of stream sediments (Lipp et al. 2020). Hence we assume a constant erosion rate ∂z/∂t = k throughout space. This simplifies Equation 1 to where |A(x)| is the total area upstream. Note that in the following sections we have dropped the x notation for simplicity. More compactly, we can write Equation 2 as D(x) = F(C).
To generate forward predictions of downstream geochemistry, the composition of sediments entering the drainage network needs to be known. To generate a continuous input field, the G-BASE samples are interpolated to a 100 × 100 m grid using continuous curvature splines (Smith and Wessel 1990). This grid is then used to parameterise the forward model. Equation 2 is implemented using LandLab and solved with the interpolated G-BASE grid in an equal-area projection.

Inverse Model
The inverse model seeks smoothly varying source region concentrations that yield low residual misfits to observed concentrations downstream. Discrete samples of composition along the river network are inverted for upstream source geochemistry, C. Following the methodology described in Lipp et al. (2021), C is discretised as a uniformly spaced (1.5 × 1.5 km) grid of vertices, which are represented as the vector C below. This grid is then upsampled using a nearest neighbour approach to the resolution of the digital elevation model (100×100 m). The upsampled grid can now be used to calculate the composition of river sediments continuously downstream using the conservative mixing model described in the previous section. To find optimal source region compositions (i.e., the smoothest best-fitting model), the misfit between observations D and theoretical downstream compositions F(C), calculated from mixing of the proposed source, are minimised. In vector-space this problem can be written as where bold letters denote vectors. This problem is likely to always be underdetermined, i.e., fewer observations than predictions (Lipp et al. 2021). Therefore, we seek the smoothest model of upstream concentrations that generate theoretical concentrations downstream that have lowest residual misfit to the CUSP observations (e.g., Parker 1994). The following objective function is minimised, where the first summand is misfit. Model roughness is calculated within the parentheses of the second summand. λ is a hyper-parameter controlling the effect roughness will have on the solution; increasing λ will create smoother solutions that fit the data less well. To calibrate λ we utilise the 'Occam's' approach to identify values that generate the smoothest models that adequately fit the data (i.e., low roughness and data misfit; see Supporting Information; Parker 1994). Formally, we identify optimal λ values as those that yield the highest curvature between roughness and data misfit.
As this inverse problem is non-linear with no obvious analytical solution, Equation 4 is minimised numerically using Powell's conjugate gradient method as implemented in SciPy (Powell 1964;Press et al. 1986). We used Imperial College's High-Performance Computing cluster to systematically explore the smoothing parameter space by running multiple inverse models concurrently. However, we note that for optimal values of λ a standard desktop computer can be used to identify best-fitting smooth solutions (i.e. run the inverse model) in less than a hour.

Quantifying excess concentrations
Calculating differences between predicted concentrations, from the forward and the inverse models, provides a means to quantify excess elemental masses that cannot be accounted for by the 'geologic' or natural contributions measured by G-BASE. The parsimonious procedure used to quantify excess elemental masses is as follows. First, the total measured suspended sedimentary flux of the Clyde river (0.11 × 10 9 kg yr −1 ; Milliman and Farnsworth 2011), is divided by the total drainage area to yield an average erosion rate of 0.059 kg m −2 yr −1 . The total eroded mass, M(x), at each point in the drainage network is calculated by integrating the average erosion rate with respect to upstream area. Finally, the following equation is used to Preprint -Quantifying excess heavy metal concentrations in drainage basins using conservative mixing models where C is the predicted fractional (i.e., 0 < C(x) < 1) concentration of a given element. This calculation is performed for both the inverse and forward model predictions. The difference between the two calculated mass fluxes gives the expected excess mass of a given element, relative to the natural baseline. Note that as the measured sedimentary flux inserted into this equation only considers the suspended load, calculated masses are an underestimate of the actual values. Alternative mass fluxes and erosion rates could be inserted into Equation 5 should they become available. Note that we are implicitly assuming that all grainsize fractions behave in the same way as the studied grainsize fraction.

Forward model predictions of river compositions:
Towards continuous natural baselines Figure 4 shows Mg, K, Sr and Mn concentrations predicted from the forward model, F(C), alongside observed concentrations, D. These elements cover a large range of concentrations, as well as different chemical affinities. K and Sr are both generally associated with felsic rocks but occur as, respectively, major and trace elements. Mg and Mn contrastingly are mafic major and trace elements. Mg, K and Sr all tightly cluster around the 1:1 line in Figure 4b, indicating minimal deviation relative to the predicted natural baseline. The root-mean-squared misfit (rms) values are low, < 0.25, for these four elements; rms, is defined as where N is the number of downstream observations. A closer inspection of the distribution of Mn misfits shows that they are roughly symmetrically distributed around a misfit of 0 (see orange histogram in Figure 4c). Figure 4d shows the predictions along the main River Clyde (lines) and the observations (symbols). These results indicate that the concentrations of these elements are largely in accordance with their predicted natural baselines. Similarly accurate predictions of As, Ca, Co, Cr, Mn, Ti and Zr concentrations are shown in Supplementary Information. Accurate predictions of downstream sediment chemistry were also found in a previous study of rivers draining the Cairngorms, Scotland (Lipp et al. 2020). Collectively they show that conservative mixing models can be used to approximate the natural baseline concentration of a wide-range of elements, should appropriate source region chemistry data be available. Observed heavy metal (Pb, Cu and Zn) concentrations, however, exceed calculated natural baselines (see Figure 5). Figure 5b shows that Cu and Zn concentrations are consistently underpredicted by 'natural' baselines calculated from the measured upstream G-BASE concentrations. Pb deviates most strongly with residuals that are broadly uniformly distributed and an rms misfit of 0.58 (see Figure 5c).

Inverse modelling to estimate excess (non-geologic) contributors to river compositions
To estimate heavy metal concentrations in the upper Clyde basin that best-fit observations downstream, we use the inverse model to 'unmix' the CUSP samples. We stress that this approach assumes conservative mixing and seeks smooth solutions (Figure 6b, e, h). Ideally, this procedure would make use of widely distributed samples to constrain the distribution of concentrations across the basin (cf. Lipp et al. 2021). As most of the samples available here lie along the main trunk of the Clyde river there is little resolving power in tributary watersheds (see Figures 16-17 in Supporting Information). We therefore focus on the best-fitting downstream river profiles for the upper Clyde. Figures 7-9 show the results for Pb, Cu, Zn.
In these examples, optimal smoothing parameter (λ) values for Pb, Zn and Cu are 10 0.2 , 10 −0.2 and 10 −0.2 , respectively (see Figure 18 in Supporting Information). Unsurprisingly, the inverse models yield lower residual misfits than the natural (G-BASE-derived) baselines. Predicted concentrations from the best-fitting, smooth, inverse models of Cu and Zn concentrations cluster around the 1:1 line (solid circles), more closely than the natural baseline (cf. empty circles; Figures 8b & 9b). The Pb residuals have a broader symmetrical distribution centred on 0 (Figure 7a).   Figure 3). Third column = average concentrations in river sediments measured in CUSP survey and their 5 th and 95 th percentile range (mg /kg; see Figures 4d & 5d). Ranges for As, Ca and Sn go to zero due to rounding of the data. Table 1 summarises measured concentrations for all elements studied. In the second column the average concentration (mg /kg) of each element in the study area measured in the G-BASE is presented with their 5 th to 95 th percentile range (see Figure  3). The third column shows average concentrations measured Preprint -Quantifying excess heavy metal concentrations in drainage basins using conservative mixing models from the 60 samples in the CUSP survey along the Clyde river and the 5 th to 95 th percentile ranges. These values provide broad constraints on the composition of river sediment. For example, deviations between G-BASE and CUSP values indicate the presence of excess concentrations of some elements downstream (e.g., Pb, Zn). We note that whilst single values such as these can be useful, it is important to recognise that river sediment geochemistry is heterogeneous and depends on the composition of upstream source areas. This heterogeneity is highlighted by the relatively large ranges of both datasets. We instead suggest that concentrations observed or calculated throughout drainage basins provide clearer guides to the provenance of material and concentrations in excess of natural background contributions.

Discussion
We have demonstrated how spot measurements of element concentrations in source regions (first order streams) can be combined with conservative mixing models to predict concentrations in rivers downstream. The main source of uncertainty for this part of the study arises from the assumption that the first-order stream samples are unaffected by human influences. This assumption is probably generally reasonable given that G-BASE sampling guidelines explicitly state that samples ought to be taken upstream of contamination sources (Johnson et al. 2005). However, it is is possible that some pollution is still present due to diffuse source (i.e., car exhausts, atmospheric fallout). Combining isotopic 'fingerprinting' with the computational methods we present could be a fruitful way to identify the origins of measured elemental concentrations. An additional source of uncertainty concerns the assumption that river sediments are well-mixed and that they are 'instantaneously' transported from upstream sources. Very recent additional elemental sources, which are yet to be well-mixed into the sediment, may be poorly characterised in the scheme we present. Whilst more complex strategies that explicitly take into account time-varying sediment transport (e.g., Palu and Julien 2019) could be considered we think that it is prudent to first test simple, conservative, mixing models. This is particularly true where detailed information about sedimentary fluxes is not available. It is encouraging that simple conservative mixing models reliably predict concentrations downstream for a suite of elements largely controlled by non-anthropogenic processes (e.g., Mg, Sr) that have diverse affinities and concentrations.
We reiterate that in the first part of this study we used G-BASE data, which seeks to avoid anthropogenic contaminants, to parameterise source region compositions. It is unsurprising therefore that elements important in industrial processes, e.g., the heavy metals Pb, Zn and Cu, have highest residual misfit to measured concentrations downstream. The concentrations of these elements tend to be under-predicted by the simple forward model of conservative mixing. We assumed that under-prediction is a consequence of additional sources of material downstream of GBASE samples that have not been accounted for. To estimate the mass and provenance of this missing material we inverted the concentrations of heavy metals from 60 localities along the Clyde river collected as part of the CUSP. By comparing concentrations calculated along the river from forward modelling of G-BASE data and from inverse modelling of CUSP data it is possible to infer sources of excess elemental masses. Figures  7d, 8d and 9d show estimated excess mass of Pb, Cu and Zn as a function of distance along the river. The inverse models predict higher concentrations of heavy metals in almost all of the source areas compared to predicted natural (geologic) baselines. An exception is for Zn at distances > 80 km, where the inverse model predicts lower concentrations than the forward model ( Figure 9c). We note that there are few observations available for inversion at distances > 95 km, as such we suggest that this result is principally a consequence of low numbers of samples. The largest increase in excess Pb occurs between ∼ 90 − 45 km. We interpret this result as an indication of a source of Pb that is not derived from natural erosion and weathering of the local bedrock (recorded by the G-BASE samples). An obvious possible source is the Lowther Hills, which have been a Pb mining area since the early Middle Ages, and have previously been identified as a candidate for Pb pollution of the Clyde (red cross in Figure 1a; Rowan et al. 1995). The Glenngonar Waters, which drains the Lowther Hills, meets the Clyde at ∼ 90 km (i.e., close to the start of the relatively rapid increase in excess Pb) and was found to have increased Pb-levels in its sediment in previous studies (Oliver and Naysmith 2011). Similar, but less pronounced steps can be seen for Cu at 45 km and 80 km for Zn.
Assuming that the sediment is well mixed with respect to upstream area (i.e., mixing is instantaneous) we estimate that 9.7 tonnes/yr (10 6 g/yr) of Pb, 1.5 tonnes/yr of Cu and 5.7 tonnes/yr are inserted into the upper Clyde sediments through anthropogenic activity. Repeat measurements to account for temporal variability in element concentrations would improve the accuracy of these estimates. The spatial resolution of the inverse model is limited by the sampling density. At present only samples from the main trunk channel were gathered. This sampling design prevents us from identifying locations of excess concentrations in most tributaries (see Supplementary Information: Figure 17). This limitation is straightforward to remedy with additional samples (see e.g., Lipp et al. 2021).
Heavy metal pollution in river sediments can be dangerous to human health and ecosystems. The "threshold effect level" (TEL) and "predicted effect level" (PEL) have been proposed as guides to ecotoxicity (see e.g., Hudson-Edwards et al. 2008;Oliver and Naysmith 2011). Concentrations of elements in river sediments below TEL are considered unlikely to pose significant hazard to aquatic organisms, while concentrations above PEL are often associated with significant adverse impacts on the ecology of rivers TEL/PEL, mg /kg:  14 consistently exceed the TEL. For Zn five samples exceed the PEL, however the inverse result never exceeds it, which we attribute to model damping. Pb consistently exceeds the PEL after the incorporation of tributaries draining the Lowther Hills. We suggest that these methodologies show promise for identifying watersheds where remediation efforts could be focused.

Conclusions
Conservative mixing models are used to calculate natural baselines for the downstream sedimentary chemistry of the Clyde river, UK. These baselines are subsequently compared to predictions from inverse modelling of measured downstream compositions to quantify heavy metal pollution in the upper reaches of the river. As far as we are aware, this approach to identifying the provenance of pollutants and determining element concentrations as continua throughout drainage basins is novel. The forward models presented use observations from first-order streams to generate continuous natural baselines. Many elements (e.g., Mg, K, Mn, Sr) with a broad range of concentrations and chemical affinities are observed to be closely aligned with their predicted baselines. Key exceptions are the heavy metals Pb, Zn and Cu, which have high misfits and are consistently underpredicted by conservative mixing models parameterised with upstream source compositions measured by the BGS in the Geochemical Baseline Survey of the Environment (G-BASE). Elemental concentrations and mixtures in source regions required to generate good fits to measured concentrations downstream along the upper Clyde river are calculated by inverting measurements made at 60 localities during the Clyde Urban Super Project (CUSP). Resultant residual misfits between observed and predicted downstream concentrations are significantly reduced. Moreover, residuals cluster around 1:1 relations if source region concentrations of Pb, Zn and Cu are considerably higher than the 'natural' concentrations measured by G-BASE. Comparing predicted heavy metal concentrations from the forward and inverse models allows us to quantify excess (non-geological) contributions to heavy metal concentrations and to tentatively infer sources of pollution. Most notably for Pb, an increase in concentration can be traced to the Lowther Hills, a historic mining region. Total pollution, i.e., cumulative masses in excess of those predicted by conservative mixing of source concentrations, can be quantified. Using this data and previously defined threshold levels, it is possible to identify foci of toxic elemental concentrations, and sources of pollutants, which we suggest could help focus remediation efforts.

Introduction
This document contains additional results from forward and inverse modelling of geochemical data. Figures 1-15 show the results from forward modelling of 15 elements recorded in the G-BASE. Figures 16-17 demonstrate the resolution of the inverse model. Figure 18 shows results from tests to identify optimal damping parameter values. . We note for completeness that 16 samples (not shown here) of grain-size fraction < 2 mm were collected from the Clyde estuary bed using Day grab, gravity corer, Mackereth piston corer or Van Veen grab. Six samples were likely also gathered on foot from the riverbank. Their compositions were measured by XRF spectrometry. The location of these samples and their compositions are given in Jones et al. (2017). Given the different sampling strategies used to collect them they are excluded from further analysis in this contribution.

Forward model predictions
In

Inverse model coverage
The area upstream of each sample site is a guide to the resolving power of the inverse model. Figure 25 shows the area between each sample and the next upstream sample (where present) or the drainage divide. The 60 CUSP samples inverted lie along the main channel, which tends to result in low model coverage (resolution) in tributary basins. Model resolution is highest in the headwaters and lowest along sparsely sampled tributaries downstream. These results are supported by 'chequerboard' tests. In these tests a known (synthetic) source is used to calculate compositions at the actual CUSP sample sites using the actual drainage network and the forward model described in the main manuscript. The source concentrations at the 60 sample sites are then inverted following the methodology described in the main manuscript. The inverse model is not preconditioned using information about the known source; in essence we have discarded the known source compositions and are attempting to recover them from only the concentrations at the 60 sample sites and the actual drainage network. The ability of the inverse model to recover source region concentrations can then be straightforwardly assessed by comparison with the known synthetic source. The results indicate that model resolution is highest in headwaters ( Figure 26). An optimised sampling strategy would include samples along the tributaries, which would drastically improve the resolution of the inverse model [cf.][their Figure 5]lipp2021. Figure 27 show an example of the relationship between model roughness and data misfit for Zn (see Equation 4 in main manuscript for details). In this series of tests, the inverse model was reran many times (with same starting conditions, input data, etc.) for a systematic sweeep of damping parameter, λ, values.
In the example shown in Figure 27, the optimal value, i.e. the one that generates the smoothest model that best fits the data is 10 −0.2 . The 'L-curves' for Pb and Cu have a similar functional form and yield optimal λ values of 10 0.2 and 10 −0.2 , respectively (see main maunscript for discussion).                 The actual drainage network was used in these calculations to test the resolving power of sample distributions and flow paths (cf. Figure 25). Note low resolution in large tributaries downstream (cf. Figure 25), and decreasing resolving power of short wavelength (< 20 km) changes in concentration. Figure 27: Trade-off between model roughness and data misfit for Zn. Example of relationship between misfit and model roughness used to calibrate λ. The optimal λ value (10 −0.2 ; smoothest model that best fits the data) is identified at the point of maximum curvature, i.e. the 'elbow' of the trade-off curve indicated by the annotated red arrow.