1 Introduction

Ultra-trace element (< 1 ng g−1) signatures in the dissolved load of river waters are valuable fingerprints of natural and anthropogenic element sources, as well as numerous soil-to-aqueous processes operating in a hydrologic catchment. Notably, the rare earth elements plus yttrium (REE + Y) (Elderfield et al. 1990; Bau and Dulski 1996; Johannesson and Hendry 2000; Lawrence et al. 2006c; Leybourne and Johannesson 2008; Kulaksız and Bau 2013; Tepe and Bau 2014; Armand et al. 2015; Duvert et al. 2015; Pédrot et al. 2015; Gill et al. 2018) and the high field strength elements (HFSE; Mo, W, Nb, Ta, Zr, Hf, Th, U) (Tosiani et al. 2004; Tepe and Bau 2014; Censi et al. 2018; Zuddas et al. 2018) can reveal insight into catchment- to global-scale chemical weathering and particle-reaction processes. These elements enter and are transported through the dissolved load in the terrestrial realm of the hydrosphere largely bound to colloids (Pokrovsky et al. 2006) with rivers ultimately constituting a major source to the oceans (van de Flierdt et al. 2004; Rickli et al. 2010). As such, the geochemical cycle and fractionation of REE + Y and HFSE (e.g. Nd vs. Hf) in the terrestrial environment exerts an important control on oceanic budgets and associated isotopic signatures available for scavenging by mineral particles or organic matter and sequestration in marine sediments (Bau and Koschinsky 2006; Godfrey et al. 2009; Rickli et al. 2009; Firdaus et al. 2011; Schmidt et al. 2014). However, the mechanisms releasing and fractionating the elements from bedrock through to riverine transport prior to the zones of estuarine mixing/submarine discharge into the oceans remain poorly understood, especially for the HFSE (e.g. Nb/Ta, Zr/Hf, Mo/W). Some documented REE + Y-HFSE fractionation patterns from crustal source to oceanic sink have been traced as far back as the Archean, such that reconstructing these cycles can reveal important aspects of the emergence and weathering of ancient crust and, consequently, the evolution of the lithosphere-atmosphere-hydrosphere through geological time (Viehmann et al. 2014, 2018).

The REE + Y and HFSE, although only sparingly soluble (Taylor and McLennan 1985), show a significantly greater magnitude of relative fractionation in aqueous systems than in solid Earth systems, due to differential mineral weathering (Patchett et al. 1984) and various ligand and electron chemistry controls (Bau 1996). By example, the isovalent and near-identical ionic radii Zr(IV)–Hf(IV) and Nb(V)–Ta(V) pairs are often colloquially referred to as “geochemical twins” due to their limited fractionation during solid Earth processes (e.g. partial melting and fractional crystallization). The Zr/Hf ratio of various compositional estimates of the upper continental crust are similar, such as MuQ (36.9) (Kamber et al. 2005), global subducted sediment (32.0) (Plank and Langmuir 1998), and the composite of Rudnick and Gao (2014) (36.4). Similarly, the Nb/Ta ratios of upper continental crust are generally between 11 and 14 (Barth et al. 2000; Rudnick and Gao 2014). By contrast, the Zr–Hf and Nb–Ta pairs show significantly fractionated and generally supercrustal ratios in the hydrosphere. The Nb/Ta and Zr/Hf ratios reported from filtered Pacific Ocean water range from 46 to 349 (mean: 184) and 13–85 (mean: 26), respectively (Firdaus et al. 2011; Niu 2012). However, aquatic Nb–Ta and Zr–Hf studies to date have focused primarily on marine environments (Firdaus et al. 2011; Schmidt et al. 2014; Firdaus et al. 2018; Censi et al. 2019). Fractionation of Zr/Hf in rivers has been more fully documented in recent years (Godfrey et al. 1996, 2008, Zuddas et al. 2017; Censi et al. 2018; Zuddas et al. 2018), but there are very limited constraints on Nb/Ta fractionation in rivers due to the scarcity of Ta data (Filella 2017). At face value, consensus points towards supercrustal ratios being generated via removal of Ta > Nb and Hf > Zr during particle interaction, but opposing trends of fractionation with changing Nb and Zr abundances points to different controlling particles, adsorption mechanisms, or biogeochemical effects that require further work to unravel (Firdaus et al. 2011; Niu 2012).

The compiled trace element chemistry of major rivers can place important baseline constraints on terrestrial fractionation processes and are used to calculate element fluxes to the oceans. However, very few rivers have been fully characterized for a comprehensive suite of dissolved ultra-trace elements (Gaillardet et al. 2014). The Ottawa River in Canada is one of the best characterized examples due to it being the source of the National Research Council – Conseil national de recherches Canada (NRC-CNRC) series of SLRS natural river water certified reference materials (CRM). The SLRS CRM are frequently used for natural water data quality assurance/quality control (QA/QC) in laboratories across the globe (Yeghicheyan et al. 2001; Barroux et al. 2006; Lawrence et al. 2006a; Dick et al. 2008; Bayon et al. 2011; Heimburger et al. 2013; Yeghicheyan et al. 2013; Hoang et al. 2019). Consequently, these data are opportunistically taken as representative of the trace element and isotopic composition of the Ottawa River for global compilation studies (Archer and Vance 2008; Vance et al. 2008; Gaillardet et al. 2014). However, the SLRS CRM series represent only a well-characterized composition of the river at a singular location and time. The only targeted trace element study of the Ottawa River with a wide spatial coverage predated the development of low-level, high-precision mass spectrometry techniques (Merritt 1975).

There are two primary facets of this study: (1) characterization of ultra-trace elements in the dissolved load of the Ottawa River and selected tributaries near to and upstream of the SLRS CRM series sampling sites; and (2) demonstration of the precision, accuracy, and portability of a method capable of rapidly determining ultra-trace elements in natural waters, including all REE + Y and several HFSE (Nb, Ta, Zr, Hf, Mo, W, Th, U). The method is applied to SLRS-6, the current SLRS CRM generation distributed by NRC-CNRC, and a suite of newly collected samples from the Ottawa River basin (ORB). This is the first study to provide baseline ultra-trace element geochemistry across the wider ORB with the aim to understand weathering source, downstream riverine fractionation, and potential roles of anthropogenic influence.

The new contributions from this study are relevant to different parts of the Earth science community. First, the Ottawa River and several of its tributaries drain Precambrian shield rocks, such that the samples capture the geochemical signature exported to rivers from the weathering of ancient evolved crust. The focus of this study is on full REE + Y patterns and anomalies (La, Ce, Eu, Gd, Y) and HFSE abundances and ratios (Nb/Ta, Zr/Hf, Mo/W, Th/U) to shed new light on the terrestrial fractionation of these elements in a silicate-dominated catchment. Notably, in addition to being useful natural process tracers, these element groups include a number of “technology-critical elements” with the potential to emerge as environmental contaminants and thus the ORB results also help improve our understanding of their environmental baselines (Filella and Rodríguez-Murillo 2017; Balaram 2019). Second, SLRS-6 is certified for only a limited number of trace elements, and certification does not include the REE + Y and most of the HFSE. Thus, the community relies on published values for inter-laboratory testing (Fisher and Kara 2016; Filella and Rodushkin 2018). Published SLRS-6 ultra-trace element data are scarce but now increasing in availability, especially for the REE (Amorim et al. 2019; Lerat-Hardy et al. 2019; Schmidt et al. 2019; Yeghicheyan et al. 2019). This study provides new SLRS-6 abundance data for 11 certified elements and 27 uncertified elements, including the HFSE that can be used as informational values. Finally, this contribution outlines a direct analysis quadrupole inductively coupled plasma mass spectrometer (Q-ICP-MS) method updated from Lawrence et al. (2006a) that can produce high-precision data for ≥ 38 elements on aliquots of water down to volumes as low as ≤ 2 mL. One of the primary barriers to understanding the REE + Y and HFSE in natural waters is their very low abundance that often requires laborious pre-concentration to overcome (Bau and Dulski 1996; Bayon et al. 2011; Viehmann et al. 2014; Fisher and Kara 2016; Viehmann et al. 2018; Hoang et al. 2019). The method outlined here requires minimal preparation and can produce data for all HFSE and REE + Y in parallel with numerous other water tracers (e.g. Sr, Ca, Mg, Na, K, Li, Rb), assuming sufficient natural abundance levels and strict control of blank through clean laboratory handling.

2 Overview and Geological Context of the Ottawa River Basin (ORB)

The Ottawa River (or “Great River” as translated from the Algonquin name Kichesippi) is classified as a lacustrine river spanning 1271 km through interconnected lakes, man-made dams/reservoirs, waterfalls, and rapids that form a significant extent of the Ontario-Québec provincial border. The Ottawa River has 28 major tributaries, 24 of which are downstream of Lake Timiskaming, and a mean discharge of 1948 m3 s−1 before draining into the St. Lawrence River (Fig. 1). The ORB covers an area of 146,334 km2 within the wider St. Lawrence river basin (Thorp et al. 2005) with a bedrock geology dominated by Archean-Proterozoic plutonic and metamorphic rocks at 88% with lesser amounts of Archean volcanic rocks (4%) and Paleozoic sedimentary rocks (8%) (Baer et al. 1978; Telmer 1997). The Precambrian bedrock that extends throughout the northern and central parts of the ORB belongs predominantly to the Archean Superior Province near its southern contact with the Proterozoic Grenville Province, whereas a smaller area of the Proterozoic Huronian Supergroup of the Southern Province is exposed in the northwest part of the catchment (Card 1990). A simplified geological map showing the extensive coverage of the Archean-Proterozoic felsic gneisses and plutonic bedrock is presented in Fig. 2 (Telmer 1997; Telmer and Veizer 1999). The upper ORB has only sparse Quaternary till cover and generally thin soils relative to the lower ORB (Shilts et al. 1987; Telmer 1997) and was also not influenced by the Champlain Sea glaciomarine sediments deposited during recession of the Laurentide Ice Sheet (Parent and Occhietti 1988; Occhietti 1989).

Fig. 1
figure 1

Simplified hydrology of the Ottawa River basin (ORB) showing the Ottawa River and main tributaries. Sampling sites along the Ottawa River transect (red squares) and 4 tributaries (Mattawa, Petawawa, Noire, Coulonge; yellow squares) for this study and for SLRS-3/-4 (near Chenaux, Ontario) and SLRS-5/-6 (near Ottawa, Ontario) are indicated. Inset map of Canada shows the position of the ORB along the Ontario-Québec provincial border. Map redrafted and modified from Telmer and Veizer (1999)

Fig. 2
figure 2

Simplified bedrock geology map of the Ottawa River basin (ORB). Map redrafted from Telmer and Veizer (1999) based on the original work of Baer et al. (1978) (bedrock geology) and Occhietti (1989) (extent of Champlain Sea incursion)

Outside of the geochemical data available for the SLRS CRM series (Archer and Vance 2008; Vance et al. 2008; Gaillardet et al. 2014), studies at the catchment scale of the Ottawa River basin to date have focused on C-O–H stable isotope biogeochemistry (Telmer 1997; Telmer and Veizer 1999, 2000). Other studies have reported a more limited major element, trace element, and Sr and S isotope dataset (Merritt 1975; Wadleigh et al. 1985; Yang et al. 1996; Telmer 1997; Rondeau et al. 2005). Data from these studies place context on the spatial distribution of natural and anthropogenic inputs to river waters of the ORB. Most important to this study is that the upper ORB has silicate rock-dominated geochemical signatures (e.g., high 87Sr/86Sr, low total dissolved solids and Ca–Mg) with limited anthropogenic influence. In the southern ORB, these signatures show a transition towards a greater influence from the Phanerozoic sedimentary rocks, agriculture, and other anthropogenic activity associated with higher population density (Wadleigh et al. 1985; Yang et al. 1996).

3 Materials and Methods

3.1 SLRS CRM and SLRS-6 Handling

This contribution considers primarily the latest four generations of the SLRS CRM, SLRS-3 to SLRS-6. All are from the Ottawa River, with SLRS-3 and SLRS-4 sampled near Chenaux, Ontario, and SLRS-5 and SLRS-6 sampled from the untreated water at the Britannia Water Purification Plant near Ottawa, Ontario (Fig. 1). In each case, river water was filtered through 0.2 µm membranes, acidified to pH 1.6 using ultrapure HNO3, stored and blended in a polyethylene tank, and subsequently bottled into individual polyethylene containers for CRM distribution.

All samples for this study were analyzed at either the Isotope Group facility at University of Tübingen (“UT Setup”) or the geochemistry facility at University of Dublin, Trinity College (“TCD Setup”). All field samples and 6 separate bottles of SLRS-6 (arbitrarily labelled UT-01, UT-02, UT-03, UT-04, UT-05, UT-06) purchased from NRC-CNRC in September 2018 were measured with the UT Setup, whereas 1 bottle of SLRS-6 (labelled TCD-01) purchased in November 2016 was measured using the TCD Setup.

For the UT Setup, a 35 mL aliquot from each SLRS-6 CRM bottle was poured from the original, well-shaken bottles into pre-leached (0.8 M HNO3 for ≥ 7 days) 50 mL polypropylene (PP) centrifuge tubes in a Class 10–100 (US FED standard class) clean laboratory. A further aliquot, referred to as UT-Mix, was prepared by combining ~ 5 mL aliquots of each individual bottle to produce a physical mixture in a separate PP test tube. The original CRM bottles were resealed for future isotopic characterization. Most measurements at UT were conducted in 2 experiments over 1 week in October 2018, with the individual bottle aliquots and the SLRS-6 Mix (k = 7) each measured in triplicate (n = 3).

For the TCD Setup, an aliquot of SLRS-6 was extracted from the original bottle prior to each experiment under a Class 10–100 hybrid fume hood in a Class 10,000 clean laboratory in a similar manner to the UT Setup. Measurements at TCD were conducted across 7 experiments from November 2016 to April 2019, one to characterize SLRS-6 comparably to the UT Setup and the remaining analyzed alongside natural water samples as part of routine QA/QC.

3.2 Field Sampling and Water Pre-treatment

All 29 field samples for this study were taken in the ORB across 10 days in late July-early August of 2018. A transect of 14 samples of the Ottawa River between Temiscaming, Ontario and Chenaux, Ontario (the ‘T-C transect’) covering ~ 300 km of flow was taken, ending at the collection site of SLRS-3/SLRS-4 (Fig. 1). These samples included both well-mixed, rapidly flowing water and gently flowing dam/reservoir water. Additional samples from 4 tributaries feeding the Ottawa River along the transect, the Rivière Coulonge (n = 2), Rivière Noire (n = 4), Petawawa River (n = 3), and Mattawa River (n = 2), as well as small lakes/ponds (n = 4) were taken for comparison. A full list of samples with their GPS coordinates is available in Supplementary Table 1. The key aspects of the chosen sample locations for this study are their position upstream from Phanerozoic sedimentary rocks of the St. Lawrence Lowlands, thicker glacial overburden with deeper soils, and regions of more extensive agricultural land use (Sect. 2). The northernmost part of the transect also starts upstream of the influence of the Champlain Sea (Fig. 2). Thus, the river water chemistry is expected to be controlled largely by the Precambrian shield rocks in the northern areas of the ORB with minimal anthropogenic influence apart from potential contributions related to current and historic mining activities.

Samples were taken from shore or boat at least 10 cm below the water surface into a collection bottle that was thoroughly rinsed by the site water. Samples were filtered through 0.45 µm nylon membranes using a ThermoFisher Nalgene system equipped with a hand or electric pump. Accordingly, the dissolved load is defined operationally in this study as the truly dissolved fraction and fine particulates/colloids passing the 0.45 µm filter. The filtration system was reused with a new filter for each sample after thorough rinsing with ultra-pure water followed by the site water. Filtered samples were transferred to low-density polyethylene (LDPE) bottles that were previously acid-leached (0.5 M HNO3 for 2 weeks), rinsed with ultra-pure water, and sealed until the point of field use. The LDPE bottles were filled with headspace only for acidification, capped, and sealed for shipping to the University of Bern. Strict cleanliness of sampling gear and careful field protocols (including pre-leaching bottles and filtration equipment and consumables, as well as a “clean hands, dirty hands” approach to field work) are vital steps in acquiring ultra-trace data from rivers, groundwater, and precipitation (e.g., Horowitz et al. 1994; Lawrence et al. 2006c; Shotyk and Krachler 2009; Shotyk et al. 2017; Gill et al. 2018). At each site, probe measurements of water temperature, pH, Eh, and dissolved oxygen were taken, as will be described and reported in a separate study. An overview of the full method from field sampling to laboratory analysis is illustrated in Fig. 3, including an example of one of the Ottawa River sampling sites (Sample RRR01), the use of a “clean hands, dirty hands” approach to water sampling and probe measurements, and the minimal laboratory steps needed to acquire ultra-trace element data.

Fig. 3
figure 3

Methodology workflow from: (1) field sampling and probe measurement applying a “clean hands/dirty hands” approach; (2) water filtration to 0.45 or 0.2 µm (example shown with a hand vacuum system and exchangeable filters) and subsequent preservation with acid, and; (3) clean laboratory preparation, ICP-MS measurement, and analysis of results

Due to international shipping logistics, samples were acidified after transfer to the clean laboratory at the University of Bern. Acidification occurred within 2 weeks of field collection using ultra-pure concentrated HNO3 to bring samples to a final acid concentration of ~ 0.5 M HNO3. No visible changes to the water colour or evidence of precipitate formation occurred prior to acidification or after acidification prior to measurement. A 10 mL aliquot of each acidified sample was transferred to a separate pre-acid leached and rinsed (as per the LDPE bottles) polypropylene (PP) centrifuge tube for shipping to the University of Tübingen for ultra-trace element analysis. All of the collected ORB samples were treated identically to the SLRS CRM samples (as per the UT Setup) from this point forward.

3.3 Laboratory Conditions and Sample Preparation for Ultra-Trace Element Analysis

Ultra-pure (≥ 18.2 MΩ) water from Millipore Milli-Q® units was used for all reagent dilution and labware rinsing in both the UT and TCD Setups. All diluted HNO3 was prepared from sub-boiling distilled concentrated stock: 3 × progressive distillation of laboratory grade acid with an Analab CleanAcids® unit in the TCD Setup and 1x distillation of pro analytical grade acid with a Savillex DST-1000® unit in the UT Setup.

For the UT Setup, a gravimetric mixture of 9.0 g of water sample and 1.0 g of a multi-element/enriched isotope (6Li, In, Re, Bi) internal standard in a 0.45 M HNO3 matrix was prepared for analysis in 15 mL PP centrifuge tubes. For the TCD Setup, a gravimetric mixture of 1.8 g of SLRS-6 and 0.2 g of a multi-element/enriched isotope (6Li, Rh, Re, Bi, 235U) internal standard in a 0.45 M HNO3 matrix was prepared for analysis in 2 mL PP micro-centrifuge tubes. The analysis solutions of both setups were thus prepared to a nominal, gravimetric dilution factor of 1.11. The different measurement volumes are based on each facility’s routine trace element workflows using alternate options on the same autosampler (Sect. 3.4) and each setup has pros and cons relative to the other. The UT Setup requires more sample volume, but offers the potential for enhanced measurement signal (e.g. increased uptake rate) and more measurement time for either better instrument counting times (increased precision) or the ability to measure more analytes. The TCD Setup consumes 5 × less sample volume, but at a sacrifice to instrument counting time and/or number of analytes that can be measured with high precision. The slightly different internal standard mixtures used at each facility is related to the desired analytes (e.g. In is measured routinely at TCD) or availability of high-purity enriched isotopes.

An additional test experiment at UT was undertaken on all 7 SLRS-6 bottles of this study (UT-01 to UT-06 and TCD-01) at a nominal, gravimetric dilution factor of 10, with 1.0 g of sample combined with 9.0 g of the multi-element internal standard prepared to the same internal standard abundances as the primary experiments.

3.4 Instrumentation

A ThermoFisher Scientific iCAP-Q ICP-MS and Elemental Scientific (ESI) SC-2 DX autosampler were used in both setups with similar sample interface configuration and data acquisition parameters (Supplementary Table 2), but with different sample introduction strategies. The UT Setup used a SC-FAST system and a 4 mL Teflon sample loop optimized for rapid, vacuum pump-driven sample uptake and washout, with the sample injected into the ICP-MS using its on-board peristaltic pump at 30 rpm. The TCD Setup sample introduction used a custom microFAST system and a 2 mL Teflon sample loop optimized for slower, syringe-controlled uptake, with the sample injected into the ICP-MS using its on-board peristaltic pump at 20–25 rpm. Both setups used either standard (STD; standard cone configuration) or high-sensitivity (STDS; high-sensitivity skimmer cone insert added) operating modes. The latter mode was only used when an appreciable boost in sensitivity from the interface region was recognized under all other instrument conditions being comparable. The use of He in the collision cell for kinetic energy discrimination (KED) was not used for this study to retain higher sensitivity for the light mass analytes (6Li, 7Li, 9Be), aiding with more accurate abundance determination of Li and Be and resulting in a more robust drift correction in the interpolated light mass region.

3.5 Measurement Method

The mass spectrometry method is built on the drift-corrected external rock calibration strategy originally developed by Eggins et al. (1997) and its adaptation for natural waters as described in Lawrence et al. (2006a). A key aspect of the method is the two-fold correction for signal drift and ionization differences between calibration standards and unknowns. An initial internal standard correction is applied to all samples relative to a pure internal standard solution measured at the onset of the experiment and using a linear interpolation for the analytes between internal standard masses. The second correction for residual drift, most prominent in the low- to mid-mass range (between 6Li and Rh or In), is applied to every analyte mass based on linear interpolation between repeat measurements of a monitor solution (bracketing every 4–6 unknowns). The experiment measure order for the new ORB samples was out of sequence from their downstream positions, and included a mix between the Ottawa River and tributary samples.

Sample washout after measurement of digested rock RM was monitored closely on all analytes to ensure return to background levels prior to injecting the natural water samples. The complete washout sequence consisted of an on-board SC-2 DX inter-sample rinse protocol with an additional 2-stage washout (blank 0.8 M HNO3 followed by 0.45 M HNO3 run as unknowns with lower acquisition time). Using this strategy, it was not necessary to run calibration standards at the end of the experiment as in Lawrence et al. (2006a). Between water samples only the on-board rinse protocol was used.

Common spectral interferences for some analyte masses, including oxides (MO+/M+), hydroxides (MOH+/M+), and isobars (M+/M+), were corrected based on routine pre-experiment measurements of production rates from pure/mixed element solutions (Dy and Ba + Nd). Other analyte interference rates were scaled linearly based on the NdO+/Gd+ production rate from these experiments to the latest bi-annual determination from pure solutions, as summarized in Supplementary Table 3 (Aries et al. 2000; Ulrich et al. 2010; Chen et al. 2017). For example, the bi-annual determination of the SmO+ on 165Ho+ contribution (as 149Sm16O+/149Sm+) alongside NdO+ on 160Gd+ (as 146Nd16O+/146Nd+) was used to scale the specific experiment’s SmO+ on Ho as: (149Sm16O+/149Sm+)experiment = (149Sm16O+/149Sm+)bi-annual × (146Nd16O+/146Nd+)experiment ÷ (146Nd16O+/146Nd+)bi-annual. Other low abundance analytes with prominent polyatomic interferences (e.g. 44Ca16O+ on 60Ni+) were either not analyzed or not reported due to the lack of He-KED use or a calculated interference correction factor. However, the method is noted to be fully customizable to run in He-KED mode to reduce polyatomic interferences. The analytes Na, Mg, K, and Ca are reported without a mathematical interference correction (e.g. 38Ar1H+ on 39K+). However, an intra-experiment blank correction for the internal standard-bearing 0.45 M HNO3 (highest for K at ~ 5–9% of the signal in waters at 1.1x dilution) accounts for most of this polyatomic background. In the case of Mg, no significant bias from an interference was suspected based on the indistinguishable abundances determined from two analyte isotopes (25Mg, 26Mg).

Calibration used the United States Geological Survey (USGS) RM W-2a (diabase) digested in HF-HNO3 and prepared to solutions with gravimetric dilution factors of 10,000–50,000 in the same carrier acid matrix and internal standard mixture as natural water unknowns. Initial preparation of the calibration standards via acid digestion is described elsewhere for the TCD Setup (Babechuk et al. 2019, and references therein) and UT Setup (Albut et al. 2018, and references therein) and in further detail in the Supplementary Materials. Calibration lines were constructed from 3 to 6 points (after full procedural blank correction) derived from measurements of variably diluted W-2a RM at the beginning of experiments. The calibration line was also cross-checked with an additional W-2a measurement at the end of an experiment, which always showed identical signal intensities to those measured at the start (after drift correction). Calibration using a single digested rock is a strategy that provides natural relative abundances of elements, forms a reasonable matrix match to river waters, minimizes washout/memory effects, and allows recalibration of results based on changing consensus of the most probable element abundances in the rock. Further, digested W-2a solutions (stored in 0.45 M HNO3) are discarded within 3–6 months and replaced with freshly prepared solutions to ensure element abundances are not modified during storage (e.g. via adsorption to PP bottle walls). Internal testing shows negligible W-2a analyte signal deterioration within ~ 6 months, likely aided by the relatively high dilution factor (1000) of the digested rock in the stock bottles, but use of the same solution for longer is not recommended. The preferred W-2a calibration abundances are provided in Supplementary Table 3, as developed, reported, and modified across previous studies (Kamber et al. 2003; Kamber 2009; Babechuk et al. 2010; Marx and Kamber 2010; Baldwin et al. 2012). The preferred W-2a W calibration abundance of 260 ng g−1 is now updated from its previous value of 240 ng g−1 based on a recent double-spike W stable isotope study (Kurzweil et al. 2018).

3.6 Estimation of Method Detection Limits and Blanks

The background equivalent concentration (BEC) is used as a proxy for the full instrument background signal. For both setups, a mean \( \left( {\bar{x}} \right) \) BEC for each element is calculated from multiple measurements (\( x_{i} \)) of an internal standard-bearing 0.45 M HNO3 carrier acid prepared alongside and identically to the samples.

$$ \bar{x} = \frac{1}{n}\mathop \sum \limits_{i = 1}^{n} x_{i} $$
(1)

For the UT Setup, a mean BEC is calculated from 40 (n) individual within-run measurements on 5 blanks over 2 experiments (1 week in 2018). The mean BEC is ≤ 1 pg g−1 for 34 elements and ≤ 100 fg g−1 for 25 of these elements. The remaining elements with a higher BEC include the major elements 44Ca (3.4 ng g−1), 25,26Mg (0.1 ng g−1), 39K (36.5 ng g−1), and 23Na (0.50 ng g−1) with higher polyatomic interference backgrounds and 86Sr (21 pg g−1) and 7Li (183 pg g−1) with a higher natural abundance or a minor impurity in the internal standard, respectively. The BEC was equivalent to ≤ 1% of the SLRS-6 signal intensity for 32 elements, between 1 and 5% for Be and Cd, between 6 and 15% for K, Sn and Ta, and 33% for Li.

For the TCD Setup, a mean BEC is calculated from 80 individual within-run measurements over 7 experiments (2016–2019) (Supplementary Table 3). The mean BEC is slightly higher for most analytes relative to the UT Setup and notably higher (> 5 ×) for a smaller group (e.g. Na, Rb Zr, Sn, Sb, Cs, Ba, La, Ce, Eu, W, Tl, Pb, U). Most of the increases probably reflect some combination of significantly poorer air purity in the mass spectrometry laboratory, greater age of the internal standards, and more diverse array of solution and laser ablation workflows on the instrument (adding more blank to the reagents, sample introduction system, etc.) in the TCD facility relative to the UT facility. The higher U BEC at TCD is specifically related to a minor impurity in the enriched 235U internal standard. Despite inter-facility differences, the mean BEC for TCD analytes is still ≤ 1 pg g−1 for 25 elements and ≤ 100 fg g−1 for 15 of these elements.

A filter for the method detection limit was taken as a signal intensity exceeding 3x the BEC value. All elements, with the exception of Li, met this requirement for SLRS-6 and most of the natural samples. However, the substantial Li background correction in both setups is considered robust due to the consistency of the 6Li internal standard impurity added equally to all samples, and thus all data for Li are reported. Blank related to handling of field samples from this study was monitored by filtering an equivalent volume of ultra-pure laboratory water through the same filter system and treating it the same way as all samples throughout the laboratory. The residual blank was typically within or near the limits of the 3 × BEC filter and generally required very minimal extra correction to data.

The primary obstacle to achieving high accuracy and precision using the analytical method applied here is minimizing the contribution of blank from sample handling (e.g. leached from filters, blank in the acid used for acidification) in addition to any blank contribution from the mixed-element internal standard and the full analytical workspace from clean laboratory to instrument space (Lawrence et al. 2006a). Additional discussion on blank levels using an earlier adaptation of this method acquiring ultra-trace REE + Y data from SLRS CRM SLRS-4 can be found in Lawrence et al. (2006b).

3.7 Estimation of Method Precision and Bias

The full method precision at both TCD and UT is routinely estimated via repeat analyses of digested, matrix-appropriate (as dictated by individual studies) rock RM and reported as the standard deviation (s) or percent relative standard deviation (%RSD) of the mean \( \left( {\bar{x}} \right) \) abundance determined for each element.

$$ s = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^{n} \left( {x_{i} - \bar{x}} \right)^{2} }}{n - 1}} $$
(2)
$$ \% {\text{RSD}} = 100 \times \frac{s}{{\bar{x}}} $$
(3)

The estimates reflect within-lab reproducibility over the period of several months to years spanning multiple studies, i.e., intermediate precision conditions (Potts 2012). Under these conditions, the %RSD on rock reference materials has been demonstrated in previous studies to be excellent (< 10% and often < 1–3%) for the analytes reported in this study in both the TCD and UT laboratories as well as others applying the same method (Kamber 2009; Marx et al. 2011; Albut et al. 2018; Rosca et al. 2018; Babechuk et al. 2019). The method precision for river water analysis under optimistic repeatability (srepeatability) conditions and more realistic intermediate (sintermediate) conditions is illustrated below with the UT (Sect. 4.1) and TCD (Sect. 4.2) SLRS-6 results, respectively. Additional discussion on the precision of natural water abundance measurements using an earlier adaptation of this method can be found in Lawrence et al. (2006a).

The method accuracy for geological materials is routinely assessed against the compiled values reported in Jochum et al. (2016) for uncertified rock RM and the certified and informational values for the rock CRM OU-6 (Potts and Kane 2005), as described recently in Albut et al. (2018). For this study, additional accuracy evaluation of rock RM is provided through UT results in the GeoPT™ proficiency testing scheme, as detailed in the Supplementary Materials. Finally, the accuracy of the SLRS-6 data in this study is assessed directly against the NRC-CNRC certified abundances, when available, or against data available in the literature for uncertified abundances (Yeghicheyan et al. 2019). In all cases, accuracy is assessed quantitatively using a %bias calculation, where \( \bar{x} \) is the data mean from this study and µ is the reference value (certified or literature).

$$ \% {\text{bias}} = 100 \times \frac{{\left(\bar{x} - \mu \right)}}{{\bar{x}}} $$
(4)

It is noted that while the accuracy of the method is anchored to the preferred W-2a calibration abundances (Supplementary Table 3) these values can be revised based on evolving consensus and previously published data can be retrospectively re-calibrated with new preferred W-2a values. The compiled W-2 data from Jochum et al. (2016) are provided in Supplementary Table 3 for comparison to the UT and TCD laboratories’ preferred values; the bias between these two data sets is < 5% for 27 trace elements, between 5 and 10% for Be, Y, Zr, Mo, and Ta, and ~ 10% for Tl and W. For W, there is less consensus on W-2a values due to fewer laboratories producing data for low W natural materials. However, excellent agreement is noted between W abundances measured with the method herein and isotope dilution analyses for several USGS RM (Babechuk et al. 2010; König et al. 2011; Babechuk et al. 2015; Albut et al. 2018; Kurzweil et al. 2018). For Tl, it is noted that the mean abundance of 536 ± 26 (2 s) ng g−1 for OU-6 recently reported in Albut et al. (2018) from this method is identical within uncertainty to the informational value of 540 ng g−1 reported in Potts and Kane (2005). Moreover, there is also good agreement between lower Tl abundance geological RM with this method and the isotope dilution data of Brett et al. (2018). Specifically, the mean (± 2 s) Tl abundances (ng g−1) of BIR-1 and BHVO-2 determined with this method were 1.27 ± 023 (Albut et al. 2018) and 1.25 ± 0.19 (Babechuk et al. 2010) and 20.17 ± 0.97 (Albut et al. 2018) and 19.33 ± 0.42 (Babechuk et al. 2010), respectively. These abundances agree within 17% of the Brett et al. (2018) values of 1.0 ± 0.4 and 23.0 ± 2.3, respectively, and overlap within uncertainty.

The results submitted to the GeoPT™ proficiency testing scheme from UT for rounds GeoPT43 (dolerite, ADS-1), GeoPT44 (calcareous shale, ShCX-1), and GeoPT45 (silicified siltstone, GONV-1) as outlined in Webb et al. (2018, 2019a and b), respectively, provide an additional accuracy assessment. These results and their implications for the method accuracy are outlined fully in the Supplementary Materials. Most importantly, data produced from all GeoPT™ materials for the primary elements of interest to this study (REE + Y, Zr, Hf, Nb, Ta, Th, U, Mo, W) show no significant bias relative to consensus values (assigned or provisional) for all results submitted and analyzed through the program (Supplementary Table 4).

4 Results

4.1 UT SLRS-6 Results, Between-Bottle SLRS-6 Consistency, and Short-Term Method Precision

The per-bottle data for the UT Setup are reported in Supplementary Table 5. Splits from the 6 separate bottles (units, k = 6) of SLRS-6 were measured in triplicate (n = 3) across 2 experiments at UT under similar instrument conditions. These experiments were designed to evaluate inter-bottle consistency using a one-way analysis of variance (ANOVA) approach (Linsinger et al. 2001), compare individual bottle means to the UT-Mix, as well as provide an estimation of precision over the duration of a short study.

An evaluation of the mean of the triplicate measurements from each of the 6 bottles (\( \bar{x} \)) indicates that there is no statistical difference with 95% confidence between the bottle means for 37 of the 38 elements. Hafnium is the only element where the null hypothesis of the ANOVA test (equivalent bottle means) is rejected, and a post hoc t test indicates that a significant difference exists only between UT-03 and UT-05. However, the low abundance of Hf and its higher blank/signal ratio suggest that apparent inter-bottle variability could also be an artefact.

A further two-tailed t-test (k = 2) comparing the mean of the triplicate measurements of the physical SLRS-6 mixture (SLRS-6 Mix; n = 3) with a mean of all measurements from individual bottles (n = 18) using a pooled standard deviation (Eq. 4) showed no significant difference at 95% confidence for 32 elements. The differences between the two means for Na, K, Sn, Pr, Ta, and Th were significant, although the percent difference between mean values did not exceed 10%.

$$ s_{\text{pooled}} = \sqrt {\frac{{\left[ {s_{1}^{2} \left( {n_{1} - 1} \right) + \ldots s_{k}^{2} \left( {n_{k} - 1} \right)} \right]}}{{\left( {n_{1} - 1} \right) + \ldots \left( {n_{k} - 1} \right)}}} $$
(5)

Overall, for all analytes measured and reported from the UT Setup in this study, there are no obvious and repeatable indications of inter-bottle heterogeneity. Given the general repeatability of the abundance measurements for the 38 elements across the 6 separate bottles and the physical mixture, a compiled full sample mean and standard deviation of all measurements (n = 21) from each bottle (k = 7) is reported in Table 1. The latter standard deviation and all others from this point forward are reported as 2 s to be at approximately 95% confidence level and consistent with the expanded uncertainty (U) with a coverage factor of 2 that is commonly used in reference material assessment studies (Linsinger et al. 2001; Heimburger et al. 2013; Yeghicheyan et al. 2013). No statistical filter for outliers is applied. The 2RSD is < 3% for 31 elements, 3–5% for 3 elements (K, Zr, Nb), and > 5% for Be (5.2%), Cd (12%), Sn (9.5%), Ta (11%), and Th (5.7%), the analytes with the lowest abundance or higher blank/signal ratios than others. The standard deviation on these data is adopted as an estimate of the method precision under repeatability conditions (srepeatability) (Potts 2012) and indicates a high level of repeatability within closely spaced (within-week) experiment batches.

Table 1 SLRS-6 element abundances (Na, Mg, K, Ca in ng g−1, all others in pg g−1) and selected anomaly and mass ratios determined with the UT and TCD Setups along with a compilation of the data from both setups

4.2 TCD SLRS-6 Results and Long-Term Method Precision

For the TCD Setup, a compilation of 42 measurements from 19 individually diluted aliquots of SLRS-6 is used to derive a mean and standard deviation. The 2RSD under these conditions is < 3% for 30 elements, 3–5% for 3 elements (Be, Nb, Tl), and > 5% for Cd (8.7%), Sn (28%), Ta (12%), and Th (12%). Of the latter elements, Cd, Sn, and Ta had higher blank/signal ratios, similar to the UT Setup data, and Th was noted in some experiments to exhibit a higher within-experiment signal drift than other heavy mass elements. These data were collected across different batches of calibration and internal standards and more variable instrument tuning and operating conditions. Thus, the standard deviation of these data provide a more realistic estimate of intermediate precision conditions (sintermediate) (Potts 2012) that are relevant to multi-purpose analytical facilities and to aqueous geochemistry studies conducted across several seasons and annual cycles. Nevertheless, these data still indicate the method offers high-precision over a longer (multi-year) measurement timeframe and illustrate that similar data quality to the UT Setup can be achieved from small (≤ 2 mL) sample volumes (Sects. 3.33.4).

4.3 Comparison of UT and TCD Results and Compiled SLRS-6 Abundance Data

The mean SLRS-6 element abundances from both facilities (n = 21, UT; n = 42, TCD) were expected to be similar due to the common analytical method (calibration, use of internal standards, external drift correction, interference correction, etc.). A comparison of the results is first presented in Table 1 through an absolute percent difference (% diff.) calculation. The percent difference between the UT and TCD mean values was < 3% for 30 elements, 3–7% for 3 elements (Na, K, Be) and > 7% for Sb (7.3%), Ta (7.7%), Tl (7.7%), Pb (11%), and Th (10%). The offsets for Na and K are likely related to small differences in background polyatomic interference production on the analyte masses (Sect. 3.5). The larger discrepancy for Be and Ta is attributed predominantly to the very low abundances of these elements (e.g. 0.06 pg g−1 for Ta) and lower and more variable signal/blank ratio. The differences between the heavy mass element (Tl, Pb, Th) abundances are less easily explained, but are suspected to be related to the use of 235U as an internal standard in the TCD Setup, which provides a more accurate deconvolution of the heavy mass instrumental drift than the projection from Bi in the UT Setup. Furthermore, in terms of Pb, there is the possibility of a difference due to variable contamination of the 2nd generation USGS geological powder RM (Woodhead and Hergt 2000; Weis et al. 2006; Kamber and Gladu 2009).

Despite the few discrepancies between the two setups, all element differences are ≤ 15% and are consistently < 10% within measurement uncertainty. Thus, both datasets are compiled into a global mean with its associated standard deviation (n = 63; Table 1). A graphical comparison of the data from each setup is illustrated in Supplementary Fig. 1, where both data sets with their respective standard deviation are normalized to the compiled mean. The compiled mean and 2 s are the recommended informational values from this study and are used from Sect. 4.5 onward unless noted otherwise.

4.4 UT Results for SLRS-6 Test Experiments at 10x Dilution Factor and Detection Limit Barriers

The SLRS-6 results determined with the UT Setup at a 10x dilution factor are reported in Supplementary Table 6. Data are reported as a mean of 3 measurements (n = 3) each of UT-01 to UT-06 and 4 measurements (n = 4) of TCD-01 in separate respective experiments undertaken ca. 6 months apart. The standard deviation (2 s, and equivalent %RSD) of the global mean of all samples (n = 22; k = 7) is used as an estimate of precision. The purpose of these data is to demonstrate the ability to retain coherent and repeatable results for most analytes at lower signal intensities (i.e. lower signal/background ratios). At this higher dilution factor, all analytes with the exception of Sn, Cd, and Ta remain above the detection limit criteria of 3x the BEC (Sect. 3.6). As expected, the precision decreases relative to primary 1.1x dilution factor analyses (Sect. 4.1). However, the precision is still <5% for 19 elements, between 5 and 10% for 11 elements, and only >10% for Be (14%), Nb (13%), Lu (11%), Hf (23%), and W (12%). Further, the mean abundances agree very well with the primary SLRS-6 compilation, with results differing by < 3% for 27 elements, between 3 and 5% for 4 elements (Na, Sb, Tl, Pb), and > 5% for K (9.5%), Be (11%), Nb (5.9%), and Th (12%), despite the decreased precision. A comparison of the compiled results at 1.1x and UT results at 10x dilution is available in Supplementary Fig. 2. Note that Li is again excluded from the BEC filter. The coherence in SLRS-6 Li abundance at 1.1x and 10x dilution (< 3% diff.), with the latter having a significantly increased relative contribution of internal standard impurity to the 7Li mass, further exemplifies the efficacy of the applied blank correction (Sect. 3.6).

The 1.1x to 10x dilution comparison illustrates that for some of the lowest abundance trace elements (e.g. the HFSE) the method will not produce data as precise or at all (below detection) for very trace element depleted waters (i.e. in some cases pre-concentration may still be necessary to acquire data for all analytes). This is evident with detection limit barriers being reached for W and Ta for some of the low-abundance tributary samples and pond/lake samples in this study. However, very low full procedural blank levels still permit quantification of a full suite of REE + Y with sufficient precision at higher dilution factors. See Lawrence et al. (2006b) for REE + Y coherency tests at even higher dilution factors.

4.5 Accuracy Comparison to SLRS-6 Certified and Literature Abundances

The SLRS-6 data reported in Table 1 include 11 (Na, Mg, K, Ca, Sr, Mo, Cd, Sb, Ba, Pb, U) of 20 certified elements and 1 (Be) of 2 elements with reference values from NRC-CNRC. This allows an accuracy assessment for 12 elements, which is done here using a %bias calculation (see Appendix A for details). Of this group, the bias is < 5% for 9 elements (Na, Mg, K, Ca, Sr, Sb, Ba, Pb, U), between 5 and 10% for Be (9.2%) and Mo (− 7.7%), and > 10% for Cd (− 12%). However, all of the elements with > 5% bias still overlap within the expanded uncertainty (or 2 s) of both data sets (Supplementary Fig. 3). Moreover, it is noted that the uncertainty on the NRC-CNRC data for Cd and Be is significant at > 20%.

The accuracy of the remaining values provided in Table 1 is assessed against the recent SLRS-6 compilation of Yeghicheyan et al. (2019), which is reported in Table 4. The Yeghicheyan et al. (2019) data were compiled from at least 15 measurements of each analyte in 1 laboratory, but for many analytes typically from ≥ 120 measurements in 9 laboratories. This %bias between the two studies is in Table 4 and a scatter plot is used to graphically compare the data (Fig. 4). The bias is < 5% for 26 elements, between 5 and 10% for Be (9.6%), Zr (5.9%), Eu (− 7.9%), and Lu (− 8.3%), and > 10% for Nb (− 68%), Cd (− 25%), Hf (− 79%), W (− 32%), Tl (− 14%), and Th (23.5%). Notably, most of the REE + Y are highly consistent between both studies (bias < 5%). The greater discrepancy for Eu is probably related to the various strategies of handling the BaO+ interference on the Eu+ analyte mass(es). The majority of analytes with a significant discrepancy are the very low abundance HFSE (Hf, Nb, Th, W). However, it is noted that the Yeghicheyan et al. (2019) HFSE data compilation is from fewer measurements and laboratories (in some cases only 1) relative to other analytes. A more detailed comparison of the HFSE and REE + Y between different data sets is considered in Sects. 4.6.14.6.3.

Fig. 4
figure 4

Comparison of compiled SLRS-6 abundance data (ng g−1 for Na, Mg, K, and Ca, and pg g−1 for all others) from this study (n = 63; Table 1) to the compilation recently published in Yeghicheyan et al. (2019) (Table 4). The full range of abundances is shown in (a) and only those elements with < 100 pg g−1 is shown in (b), as highlighted by the shading in (a). The solid and stippled lines represent equal abundances and ± 10% difference, respectively. Error bars are the standard deviation (2 s) or expanded uncertainty (U) for the respective datasets. In (b) the elements with a bias greater than ± 10% between each study are identified

4.6 Ottawa River Basin Field Samples and Comparison to the SLRS CRM Series

The ultra-trace element data for the Ottawa River samples, selected tributaries, and ponds/lakes are reported in Table 2, Table 3, and Supplementary Table 7, respectively. The pond/lake data are secondary to the study and not considered further in the results and discussion. All elements measured in the 14 new Ottawa River samples and the 4 tributaries (Mattawa, Petawawa, Noire, Coulonge) are compared to SLRS-6 in terms of relative abundance (percent difference) in Fig. 5. The Ottawa River samples are displayed using a box and whisker plot (Fig. 5a) that provides a visual illustration of downstream element behaviour. Tight data dispersion (e.g. K, Sr, Ba) elucidate elements with conservative downstream behaviour, whereas the widest dispersion elucidates elements with the highest downstream modification (e.g. Zr, Nb, Ta, Hf, Th). Within the 14 Ottawa River samples, 1 (KLR01) is significantly higher in most trace elements. However, the median and mean abundances of most trace elements are within ± 25% of SLRS-6 with the notable exceptions of Sr, Sn, Sb, Ba, and Pb, with mean relative abundances (± 2 s) of − 39 ± 9%, 849 ± 1705%, − 81 ± 5%, − 31 ± 15%, and − 67 ± 28%. The higher Sr and Ba in SLRS-6 relative to the upstream samples of this study probably reflects the additional source of these elements from Paleozoic sedimentary rocks in the southern ORB (Figs. 1 and 2). This is consistent with the higher major element abundances in SLRS-6 (chiefly Na, Ca, and Mg) and previous observations of increasing carbonate weathering flux in the southern ORB (Telmer 1997; Telmer and Veizer 1999). Similarly, the higher Sb and Pb in SLRS-6 could point to increasing anthropogenic inputs in the more populated areas of the southern ORB. Alternatively, some of these elements could have been added during handling of the SRM before distribution. Tin has significantly higher abundance and variability in the new Ottawa River samples relative to SLRS-6. The latter trends appear to reflect some combination of the low natural Sn abundance, higher uncertainty on the analytical measurements in SLRS-6 (23%), and potentially other analytical biases, such as filtration and handling blank differences or undetected bias from analytical interferences. The mean element abundances of each tributary show an expectedly greater difference relative to SLRS-6 (Fig. 5b), which includes, in general, lower trace element abundances, apart from the REE + Y.

Table 2 Ottawa River element abundances (Na, Mg, K, Ca in ng g−1, all others in pg g−1) along with selected anomaly and mass ratios
Table 3 ORB tributary element abundances (Na, Mg, K, Ca in ng g−1, all others in pg g−1) along with selected anomaly and mass ratios
Fig. 5
figure 5

Comparison of new ORB water data to SLRS-6 via percent difference. The data are divided into the Ottawa River T-C transect (a) and Ottawa River basin tributaries (b). The Ottawa River sample data (n = 14) are presented in box-whiskers, with boxes divided at the median and extending to upper and lower quartiles and whiskers extending to the maximum and minimum data. The mean is shown as a star in each box. The tributaries are shown with only the mean percent difference for each tributary

Of the full element list produced for SLRS-6 and the ORB samples, the remaining results and discussion focus on the REE + Y and HFSE (Zr-Hf-Nb–Ta–Th-U–Mo-W). Abundances are all reported as a mass fraction (in pg g−1) and element ratios are mass/mass. Mean element ratios are calculated from individual measurements when possible or otherwise calculated from the reported mean abundance in a literature source. All uncertainties on calculated means are expressed as 2 s.

4.6.1 REE + Y

The total abundance of the REE (ΣREE) is expressed as the sum of all lanthanide abundances from La-Lu. The REE + Y data of waters are typically evaluated normalized to an upper continental crust composite, providing a baseline to evaluate bulk REE + Y fractionation in terms of pattern enrichment or depletion relative to the crust and determine anomalous element behaviour. Normalization is performed here using the alluvial sediment composite ‘Mud from Queensland’ (MuQ) (Kamber et al. 2005). Normalized ratios are reported as REEn/REEn and REE anomalies (Lan/Lan*, Cen/Cen*, Eun/Eun*, Gdn/Gdn*) are defined as departures from a smooth inter-element pattern on MuQ-normalized log-linear REE + Y plots. The REEn* values are calculated using a geometric mean, where Lan* = Prn x (Prn/Ndn)2, Cen* = Prn x (Prn/Ndn), Eun* = (Sm2n x Tbn)1/3, Gdn* = (Tb2n x Smn)1/3 (Lawrence et al. 2006c; Lawrence and Kamber 2006). The Y anomaly is instead considered in terms of the Y/Ho mass ratio.

Normalized REE + Y plots of SLRS-6 data from this study, Schmidt et al. (2019), and Yeghicheyan et al. (2019) are reported in Fig. 6a. The REE + Y patterns of all data sets overlap within uncertainty, but show slightly more deviation in the HREE relative the LREE. It is noted that the ‘smoothness’ of the non-anomalous (or potentially non-anomalous) elements in the normalized pattern (i.e. excluding La, Ce, Eu, Gd, Y, Lu) also attests to the quality of the data and accuracy of inter-element REE ratios (Lawrence and Kamber 2006). A comparison of SLRS-6 REE + Y patterns to those of SLRS-5 (Yeghicheyan et al. 2013) and SLRS-4 (Lawrence et al. 2006a) is presented in Fig. 6b. Despite changes in the ΣREE between each of these SLRS CRM generations, with SLRS-4 > SLRS-6 > SLRS-5, the shape of the normalized REE + Y patterns are nearly parallel. The inter-generation pattern similarity is also evident through Prn/Ybn ratios, with SLRS-6 (2.09 ± 0.04) being similar to those of SLRS-5 (1.86) and SLRS-4 (2.36). All three generations also share similar REE + Y anomalies. For example, the Cen/Cen* value of SLRS-6 from this study (0.575 ± 0.011) is identical to that from Yeghicheyan et al. (2019) and matches well with the calculated values of 0.62 and 0.59 from SLRS-5 and SLRS-4, respectively. Similarly, the Y/Ho ratio from this study (30.29 ± 0.76) is similar to the ratio of the SLRS-6 (29.8), SLRS-5 (30.3), and only slightly higher than SLRS-4 (28.0).

Fig. 6
figure 6

MuQ-normalized REE + Y patterns for SLRS-6 (a), selected data for SLRS-4 and SLRS-5 for comparison to SLRS-6 (b), new samples from the Ottawa River from this study (c), and the same data from (c) scaled to the ΣREE of SLRS-6 to illustrate pattern similarity apart from variations in Ce, Eu, and Y (d). The SLRS-6 data from this study (Table 1) are compared to those of Yeghicheyan et al. (2019) (Table 4) and Schmidt et al. (2019) (excluding Y) with respective uncertainties. Data for SLRS-4 (filled diamonds) and SLRS-5 (filled hexagons) are from Lawrence et al. (2006a) and Yeghicheyan et al. (2013), respectively

One exception to REE + Y similarity between the latest 3 SLRS generations is a positive Sm anomaly observed only in SLRS-4 (Fig. 5b), as previously documented across several independent studies (Yeghicheyan et al. 2001; Lawrence et al. 2006a; Lawrence and Kamber 2007). Yeghicheyan et al. (2013) attributed this Sm over-abundance to uncorrected interferences on the Sm analyte mass. However, there are not many documented polyatomic or doubly-charged interferences on the most common Sm analyte masses (149 in this study). Lawrence and Kamber (2007) instead advocated for a contamination origin, noting that other NRC-CNRC natural water CRM (i.e. SLEW, CASS, NASS) contained a similar Sm over-abundance that was reproducible in different laboratories and across different generations of these CRM. Thus, although anthropogenic Sm anomalies linked to industrial waste have been documented in river waters (Kulaksız and Bau 2013), the lack of a Sm anomaly in the later generation SLRS-5 and SLRS-6 CRM appears to converge on Sm having been contaminated only for SLRS-4 in this CRM series, at some point during or after field sampling and before bottling and distribution. The new ORB samples from this study are also devoid of Sm anomalies, yet show highly similar REE + Y patterns to the SLRS CRM (Figs. 6 and 7; see below). Collectively, these new Ottawa River observations strengthen the case for a highly selective Sm contamination that is restricted to older generation NRC-CNRC natural water CRM.

Fig. 7
figure 7

MuQ-normalized REE + Y pattern for the Rivière Noire (a), Rivière Coulonge (b), Petawawa River (c), and Mattawa River (d) tributaries in the ORB. The pattern of SLRS-6 (pink line; Table 1) is shown in each panel for comparison

The ΣREE of the 14 T-C transect samples taken from the Ottawa River (921–1610; mean: 1156 ± 393) scatter closely to that of SLRS-6 (ΣREE of 973 ± 5), although there is a broadly defined downstream decrease in abundance (Fig. 8). All of the T-C transect samples show highly similar REE + Y patterns to the SLRS CRM (Fig. 6c), which is expressed visually by mathematically scaling the T-C transect sample ΣREE equivalent to be equivalent to that of SLRS-6 (Fig. 6d). Further, the Ottawa River samples have similar Prn/Ybn ratios, with T-C samples ranging from 1.95 to 2.23 (mean: 2.11 ± 0.17) and overlapping that of SLRS-6 (2.09 ± 0.04). Inter-sample REE + Y variability is largely restricted to Eu and Ce anomalies, and to a lesser extent, Y anomalies (Figs. 6d and 8). Both Eun/Eun* and Cen/Cen* values have true anomalies that decrease gradually downstream along the T-C transect with a range of 0.974–0.853 and 0.715–0.545, respectively, that end near to the SLRS-6 values (Eun/Eun*: 0.854 ± 0.022; Cen/Cen*: 0.575 ± 0.011). The Lan/Lan* values (0.991–1.034; mean: 1.015 ± 0.024) along the T-C transect scatter close to the SLRS-6 value (1.024 ± 0.033), but do not reveal the development of an anomaly or a downstream trend. The Gdn/Gdn* values (1.050–1.098; mean: 1.072 ± 0.027) scatter close to the SLRS-6 value (1.073 ± 0.028) with all revealing small positive anomaly development, but no defined downstream trend. The Y/Ho ratios along the T-C transect (28.26–29.79; mean: 29.07 ± 0.76) reveal a small positive Y anomaly without any downstream trends, but are slightly lower than in SLRS-6 (30.29 ± 0.76). The 4 tributaries have more variable REE abundances, normalized REE + Y pattern slopes, and REE anomalies relative to each other and the Ottawa River samples (Fig. 7), including reaching lower Ce anomalies (Petawawa Cen/Cen*: 0.443–0.477), lower Eu anomalies (Petawawa Eun/Eun*: 0.720–0.733), more positive Gd anomalies (all tributaries Gdn/Gdn*: 1.08–1.12) and steeper, more LREE-enriched normalized REE + Y patterns (Colounge Prn/Ybn: 2.867–2.875). With the exception of the 2 Mattawa River samples, there is generally a high degree of consistency in REE characteristics within each tributary (Fig. 8).

Fig. 8
figure 8

Spatial changes in ORB REE + Y characteristics: a ΣREE; b Prn/Ybn; c Lan/Lan*; d Cen/Cen*; e Eun/Eun*; f Gdn/Gdn*; g Y/Ho. Samples are divided into SLRS-6 (Table 1), new Ottawa River samples (Table 2), and tributaries (Table 3). Data in each plot are organized from upstream (top) to downstream (bottom) for each respective river

4.6.2 Zr-Hf and Nb–Ta

The Zr abundance of SLRS-6 from this study (65.8 ± 1.6) overlaps with the abundance from Yeghicheyan et al. (2019) (62 ± 11). From previous SLRS generations, reported abundances of Zr were 30 ± 7 (Yeghicheyan et al. 2019), 20 ± 30 (Yeghicheyan et al. 2013), and 35.7 ± 3.4 (Hoang et al. 2019) for SLRS-5 and 61.8 in SLRS-3 (Firdaus et al. 2008). The Hf abundance of SLRS-6 from this study (2.00 ± 0.06) is significantly lower than the abundance from Yeghicheyan et al. (2019) (9.5 ± 0.4). From previous SLRS generations, reported abundances of Hf were 1.12 ± 0.20 in SLRS-5 (Hoang et al. 2019) and 0.87 in SLRS-3 (Firdaus et al. 2008). The Zr/Hf ratio of SLRS-6 in this study (32.84 ± 1.07) differs substantially from that of Yeghicheyan et al. (2019) (6.5), but overlaps within uncertainty of the ratio for SLRS-5 (31.88) from Hoang et al. (2019) despite the lower SLRS-5 Zr-Hf abundances. The Zr/Hf ratio of SLRS-3 from Firdaus et al. (2008) is higher at ~ 71. Notably, the Zr/Hf ratios reported for SLRS-5 and SLRS-6 are within the range expected for upper crustal rocks and were reproducible in the newly collected Ottawa River samples (see below). Collectively, these observations suggest that the lower Hf abundance for SLRS-6 reported here is more accurate than that of Yeghicheyan et al. (2019), but further independent tests are needed. The new Ottawa River samples of the T-C transect have a relatively consistent Zr/Hf ratio (32.46 ± 5.07) that appears to show some deviation towards lower ratios in the samples furthest downstream (SBB01, PDF01 with 26.86–28.74) (Fig. 9). Across the transect, the Zr and Hf abundances show a tightly coupled downstream depletion. The Zr and Hf abundances of the 4 tributaries are lower than in the Ottawa River samples, and Zr/Hf ratios are relatively consistent within each tributary and between each other, but significantly lower than the Ottawa River with a mean of all 11 samples at 19.19 ± 3.21 (Fig. 9).

Fig. 9
figure 9

Spatial changes in ORB HFSE ratios (a: Zr/Hf; b: Nb/Ta; c: Mo/W; d: Th/U) divided into SLRS-6 (Table 1), new Ottawa River samples (Table 2), and tributaries (Table 3). Data in each plot are organized from upstream (top) to downstream (bottom). Shaded horizontal areas represent upper continental crust Zr/Hf and Nb/Ta ratios

The Nb abundance of SLRS-6 from this study (2.60 ± 0.10) is significantly lower than that reported by Yeghicheyan et al. (2019) (8.1 ± 5.7), albeit still within their uncertainty. From previous SLRS generations, reported abundances of Nb were 6.5 ± 4.8 (Yeghicheyan et al. 2019), 3.6 ± 1.6 (Yeghicheyan et al. 2013), 3.8 ± 0.6 (Heimburger et al. 2013), 3.70 ± 0.20 (Hoang et al. 2019), and 2.8 (Filella et al. 2014) for SLRS-5, 3.89 ± 0.26 (Filella and Rodushkin 2018) for SLRS-4, and 2.9 (Firdaus et al. 2008) for SLRS-3. There are no other SLRS-6 Ta data available at present for a comparison to the abundance of 0.060 ± 0.008 reported here. From previous SLRS generations, reported abundances of Ta were 0.062 ± 0.010 (Hoang et al. 2019) in SLRS-5, 0.15 ± 0.03 (Filella and Rodushkin 2018) in SLRS-4, and 0.25 (Firdaus et al. 2008) in SLRS-3. The SLRS-6 Nb/Ta of 43.4 ± 5.5 from this study is lower than in SLRS-5 (59.7), but higher than in SLRS-4 (26.4 ± 5.9) or SLRS-3 (11.4), based on the aforementioned studies. The new Ottawa River samples of the T-C transect show a spread in Nb/Ta ratios from 20.0 to 30.5 with a gradual downstream increase (Fig. 9). Across the transect, Nb and Ta abundances decrease progressively downstream similar to Zr-Hf. The Nb and Ta abundances of the 4 tributaries are lower than in the Ottawa River samples and are generally shifted towards higher Nb/Ta ratios (29.3–45.2) that scatter near the value in SLRS-6 (Fig. 9). The lower Ta abundances of the tributaries were closer to the detection limit than the Ottawa River, and in the case of the Petawawa River and ponds/lakes (Supplementary Table 7) the samples were below the detection limit filter. The determination of Ta in natural waters is notoriously challenging even with carefully controlled blanks and care to minimize instrumental memory effects (Filella and Rodushkin 2018). Further work using different analytical strategies to better constrain Ta (and other low-solubility, low-abundance HFSE) in SLRS-6 and other natural water CRM is clearly needed. Nevertheless, the method outlined herein has excellent Ta accuracy for rock RM, and a good repeatability of SLRS-6 Ta abundances and Nb/Ta ratios was found across two separate facilities (Table 1). The new Ottawa River samples have Ta abundances at or above that in SLRS-6 such that the Nb–Ta uncertainty for SLRS-6 should provide a reasonable estimate for these samples. This is not the case for the tributary samples, where precision is poorer at lower analyte levels and the lower signal/blank ratios can account for more of the Nb/Ta scatter.

4.6.3 Mo-W and Th-U

The Mo abundance of SLRS-6 from this study (198.4 ± 5.5) overlaps with the abundance of Yeghicheyan et al. (2019) (196 ± 18). From previous SLRS generations, reported abundances of Mo were 209 ± 15 (Yeghicheyan et al. 2019), 220 ± 10 (Yeghicheyan et al. 2013), and 227 ± 8 (Hoang et al. 2019) for SLRS-5 and 201 (Firdaus et al. 2008) for SLRS-3. Notably, these measurements all closely approach the certified NRC-CNRC Mo abundances for the different generations. The W abundance of SLRS-6 from this study (11.24 ± 0.25) is significantly lower than in Yeghicheyan et al. (2019) (16.5 ± 0.4). From previous SLRS generations, reported abundances of W were 10 ± 3 (Yeghicheyan et al. 2019) and 14 ± 18 (Yeghicheyan et al. 2013) for SLRS-5 and 4 (Firdaus et al. 2008) for SLRS-3. The Mo/W ratio in SLRS-6 determined in this study (17.65 ± 0.62) is higher than the ratio from Yeghicheyan et al. (2019) (12). The SLRS-5 Mo/W ratios are only slightly higher at 20–21 and the SLRS-3 ratio is higher at 50, based on the aforementioned studies. The new Ottawa River samples of the T-C transect show a range of Mo/W ratios (11.97–23.46; mean: 17.77 ± 7.50) that scatter around the SLRS-6 value. Spatially, the Mo/W decreases gradually downstream from the uppermost Timiscaming samples (TMU01, TMD01, TML01) to the 3 samples near Rolphton (RJR01, RJR02, and RJL01) at which point the remaining downstream samples switch to gradually increasing Mo/W. The tributary samples all have lower Mo-W abundances and higher Mo/W ratios (ranging from 43.63 to 74.52) relative to the Ottawa River samples.

The Th abundance of SLRS-6 from this study (19.8 ± 2.7) overlaps within uncertainty with the abundance of Yeghicheyan et al. (2019) (16 ± 7). From previous SLRS generations, reported abundances of Th were 4.5 ± 5.4 (Yeghicheyan et al. 2019), 13.6 ± 3.3 (Yeghicheyan et al. 2013), and 14.6 ± 0.4 (Hoang et al. 2019) for SLRS-5. The U abundance of SLRS-6 from this study (70.4 ± 1.4) overlaps within uncertainty with the abundance of Yeghicheyan et al. (2019) (67 ± 3). From previous SLRS generations, reported abundances of U were 91 ± 12 (Yeghicheyan et al. 2019), 93 ± 15 (Yeghicheyan et al. 2013), and 93.9 ± 2.2 (Hoang et al. 2019) for SLRS-5. Notably, these measurements all closely approach the certified NRC-CNRC U abundances for the 2 SLRS generations. The Th/U ratio for SLRS-6 from this study (0.28 ± 0.04) is similar to the ratio from Yeghicheyan et al. (2019) (0.23 ± 0.16), but both are higher than the range of ratios for SLRS-5 (0.05–0.16) based on the aforementioned studies. The new Ottawa River samples of the T-C transect have Th/U ratios similar to SLRS-6, ranging from 0.17 to 0.37 (mean: 0.28 ± 0.15), with a broad downstream decrease (Fig. 9). The Th/U ratios of the Mattawa River (0.41–0.22) and Petawawa River (0.30–0.41) tributaries are near to or slightly higher than the Ottawa River, whereas those of the Rivière Noire (0.74–0.94) and Rivière Coulonge (0.77–0.85) tributaries are significantly higher.

5 Discussion

5.1 Filtration, Preservation, and Storage Effects on Sample Chemistry

The ORB samples from this study were filtered to 0.45 µm compared to the 0.2 µm used for the NRC-CNRC SLRS series CRM. The filtration difference could produce an abundance bias through colloids or particulates between the filter membrane sizes; however, a high degree of consistency in the chemistry between the new Ottawa River samples and SLRS-6 was observed (Sect. 4.6, Fig. 5). Importantly, this includes no major offsets in the relative REE + Y patterns (Fig. 6), consistent with the findings of Lawrence et al. (2006b), or in HFSE ratios (Nb/Ta, Zr/Hf, Nb/Ta, Th/U, Mo/W). Thus, even with some abundance bias possible (most prominently for the low-solubility Nb–Ta-Zr-Hf–Th with a greater affinity for colloids), inter-element ratios appear coherent and, accordingly, any influences from the filtration difference are considered negligible relative to other analytical effects or natural processes in the catchment. Modification of the new ORB sample chemistry after filtration and before acidification and measurement can also be ruled out on similar grounds, and is constent with the lack of visual change to the sample solutions during this timeframe. Conversely, the consistency observed points to minimal contamination of the studied elements during the handling of the SLRS-6 CRM before distribution (note the potential exception of Pb, Ba, and Sb; Sect. 4.6).

Several laboratories retain bottles of SLRS CRM for QA/QC beyond their official shelf-life period (SLRS-6 due to expire in September 2020). The close match in chemistry between SLRS-6 and the new samples from this study also demonstrates that the REE + Y and HFSE signatures appear to be retained in the CRM containers throughout the existing SLRS-6 shelf-life (ca. 4–5 years). Further, no gradual signal loss was noted across the longer measurement timeframe (~ 30 months) of the single SLRS-6 bottle (UT-01) using the TCD Setup. The modification of sparingly soluble HFSE (Zr, Hf, Nb, Ta, Th) abundances in nitric acid solutions over longer time periods is possible, but there are insufficient data to quantify this at present. One potential exception is Th. The SLRS-5 Th abundances and Th/U ratios reported in Yeghicheyan et al. (2019) are lower and less precise in comparison to the earlier Yeghicheyan et al. (2013) study, yet U abundances were similar. These observations suggest variable removal of Th but not U from the SLRS-5 CRM solution (e.g. via adsorption to container walls) during the longer aging period of the bottles between the latter studies. It is also noted that the high variability of Ta across different SLRS CRM series could encompass similar Ta removal in distributed bottles over time, but this can’t be separated from analytical or natural effects with existing information.

5.2 REE + Y Pattern Implications for Element Sourcing, Riverine Fractionation, and Long-Term Catchment Flux Consistency

The “dissolved” REE + Y are affiliated primarily with the colloidal fraction of river and ground waters (Ingri et al. 2000; Andersson et al. 2006; Pourret et al. 2009; Pokrovsky et al. 2010) and normalized REE + Y patterns of this load are variably fractionated relative to presumed rock sources and bed and suspended loads (Sholkovitz 1995). Following release into the hydrosphere, aqueous REE + Y geochemistry is controlled by REE-complexation and the extent of colloid interaction with particulates, which is in turn driven by numerous physicochemical parameters such as pH, redox state, and dissolved organic matter content. As such, pH and ionic strength are informative first proxies for monitoring conditions that favour particulate coagulation and enhanced REE + Y removal (Elderfield et al. 1990; Sholkovitz 1995; Johannesson et al. 2006; Lawrence and Kamber 2006; Pourret and Tuduri 2017) or the pH-dependent stability of REE-complexes (Millero 1992; Johannesson et al. 1995; Quinn et al. 2006; Noack et al. 2014). Dissolved REE abundances and REE + Y patterns in an individual catchment are also amenable to seasonal changes inducing variable runoff and water–rock interaction and ligand chemistry (Dia et al. 2000; Hagedorn et al. 2011; Gill et al. 2018).

Across the Ottawa River T-C transect, the full range in ΣREE is equivalent to a 43% change in concentration, indicating downstream REE + Y removal processes and/or dilution effects. However, no clear correlation of the Ottawa River ΣREE with field parameters (e.g. pH, Eh) is evident, making it difficult to elucidate colloid-particulate controls relative to dilution and mixing effects from tributaries. Nevertheless, all samples show a coherent REE + Y pattern (excluding minor changes in Ce and Eu) indicating that ΣREE changes occur with only minor inter-element fractionation, echoing earlier studies showing downstream REE + Y pattern consistency (Lawrence et al. 2006b). Short-term seasonal changes in REE chemistry of the ORB is untestable in this study due to the limited sampling period. However, the new Ottawa River observations coupled with the SLRS CRM series data illustrate that the primary REE + Y pattern characteristics of the Ottawa River (LREE > HREE enrichment, deep negative Ce anomaly, minor negative Eu anomaly, minor positive Gd and Y anomalies) are retained over a decadal timescale (SLRS-3/-4 CRM sampled in the 1990s). Further, the similar observations from 4 tributaries of the Ottawa River demonstrate that these REE + Y characteristics are a wider feature of the ORB with only minor localized changes in REE + Y patterns (e.g. slightly greater LREE/HREE enrichment, magnitude of Ce anomaly). Remaining focus is on outlining potential weathering, source rock, and riverine effects on the main features of the primary normalized REE + Y pattern in the ORB.

The minor positive Gd and Y anomalies in the Ottawa River are generated upstream of the central ORB and appear to be natural features. Gadolinium anomalies have been linked to hospital effluents (e.g. Bau and Dulski 1996; Elbaz-Poulichet et al. 2002; Kulaksız and Bau 2013; Lerat-Hardy et al. 2019), but the minor positive anomalies of the ORB are consistent throughout a wide stretch of the Ottawa River (from the uppermost T-C transect to the sampling point of SLRS-6 outside of Ottawa) and present in tributaries that are not downstream of any anthropogenic sources. The ORB Gdn/Gdn* values are also within the range of those from streams in South East Queensland that were interpreted to be natural positive anomalies by Lawrence et al. (2006b). Both of these catchments apparently have processes in operation favouring preferential Gd release and/or complexed aqueous stability, but which remain poorly understood. Similarly, the minor Y anomalies cannot be readily attributed to a source of marine carbonate (Ryan et al. 2018), with only minor inliers of Paleozoic carbonate present upstream, or marine phosphate-bearing agricultural fertilizer (Lawrence et al. 2006b; Marx et al. 2010; Kechiched et al. 2020). A potential exception may be revealed through the chemistry of the SLRS-5 and SLRS-6 CRM, which were sampled downstream of the aforementioned sources and have slightly higher Y/Ho ratios than SLRS-3/-4 and the new upstream ORB samples. The minor positive Y anomalies throughout the remaining ORB are also in line with observations from other rivers (Lawrence et al. 2006b; Leybourne and Johannesson 2008). In these cases, the anomalies are likely inherited from chemical weathering processes, where the release of fluids from weathering profiles with Y/Ho ratios greater than the source rock has been attributed to preferential Ho scavenging by (oxy)(hydr)oxides (Thompson et al. 2013; Babechuk et al. 2015), similar to observations from marine environments (Bau 1999; Bau and Koschinsky 2009).

Riverine Ce anomalies are widely accepted to record oxidative processes in soils and shallow groundwater that involve oxidation of Ce(III) to Ce(IV) and its preferential scavenging on soil or bedrock fracture Fe–Mn (oxy)(hydr)oxides (De Carlo et al. 1997; Laveuf and Cornu 2009; Pédrot et al. 2015; Yu et al. 2017). The Cen/Cen*values in the ORB rivers (0.443–0.715; mean of 0.567 ± 0.135) are very low and present in the upstream samples, indicating that preferential solid-phase Ce(IV) > REE(III) retention begins early in the catchment REE + Y cycle. Efficient oxidative Ce(IV) scavenging can be inferred during weathering, possibly amplified by repeated Fe–Mn (oxy)(hydr)oxide precipitation-dissolution cycles (Yu et al. 2017), a greater controlling role of (oxy)(hydr)oxides relative to soil organic matter (Pédrot et al. 2015), or the chemical resistance of Ce(IV)-bearing minerals (e.g. zircon) relative to other REE-bearing accessory minerals (e.g. apatite) in felsic bedrock (Fu et al. 2019). Independent of the initially negative Ce anomaly, the downstream decrease in Cen/Cen* along the Ottawa River transect and variable Cen/Cen* in tributaries point to riverine processes modifying Ce anomalies. Strong positive correlations of Cen/Cen* with HFSE (i.e. Zr, Hf, Nb, Ta, Th) abundances are apparent across all ORB samples, suggesting either progressive mixing with tributary waters characterized by lower HFSE/more negative Ce anomalies downstream in the Ottawa River or that the removal of colloidal HFSE and REE + Y favours Ce(IV) > REE(III). Similar riverine effects may also account for the increasingly negative Eun/Eun* moving downstream.

Streams inherit the REE + Y signatures of rocks in their host catchment to varying extents even amidst the complex processes that can shape freshwater REE + Y patterns along the continuum from weathering profile to water (Biddau et al. 2002; Lawrence et al. 2006c; Sklyarova et al. 2017). In this regard, the extreme LREE > HREE enriched normalized REE + Y patterns from the ORB (mean Prn/Ybn of 2.23 ± 0.62) are atypical relative to most river waters. The normalized REE + Y patterns of most major rivers are flatter to slightly LREE-depleted, reflecting a closer geochemical similarity to the average upper continental crust and the tendency for LREE > HREE removal from the dissolved load (Elderfield et al. 1990). To illustrate this, the ORB Prn/Ybn ratios are compared to the world river compilation of Gaillardet et al. (2014) (Fig. 10) with a mean of 1.11 ± 1.24 (n = 17, from samples with both Pr and Yb data excluding the Ottawa River). Accordingly, the ORB river waters appear to inherit a disproportionate supply of REE + Y from the highly LREE > HREE enriched felsic plutonic and metamorphic rocks in the catchment (Feng and Kerrich 1992). The REE + Y supply to ORB waters could be driven by the preferential dissolution of LREE-enriched accessory minerals (Harlavan and Erel 2002). However, the chemical weathering of coarse-grained intermediate-felsic rocks tends to result in the accumulation of LREE > HREE in weathering residues even after considering the complex controls from protolith accessory mineralogy, adsorption capacity of pedogenic mineralogy, and profile maturity (Nesbitt and Markovics 1997; Aubert et al. 2001; Yusoff et al. 2013; Fu et al. 2019). Thus, the geologic imprint from the bedrock likely overwhelms weathering or riverine effects on the REE + Y patterns in the ORB waters. Additional support for this comes from testing normalization of the waters to potential LREE > HREE-enriched source rocks in the upper ORB (Feng and Kerrich 1992), which flattens the REE + Y patterns (Biddau et al. 2002). The imprint of the silicate rock-dominated REE + Y signature in the ORB is retained downstream of the Champlain Sea sediments and other Paleozoic sedimentary rocks (i.e. southern most samples of this study and SLRS-5/-6 sampling area).

Fig. 10
figure 10

Histogram of Prn/Ybn ratios (MuQ-normalized) from the world river compilation of Gaillardet et al. (2014) and the ORB (including Ottawa River and tributaries) to illustrate the atypical LREE > HREE enrichment in the latter

The propensity for geologic signatures to be retained across a major river catchment is relevant to deciphering variations in different shallow marine waters inferred to have a strong continental REE + Y source inheritance (e.g. Molina-Kescher et al. 2018). Similarly, these observations are relevant to REE + Y hydrosphere cycling on the Archean-Proterozoic Earth surface where a dominating continental solute budget is recorded in near-shore sedimentary deposits (Alexander et al. 2008). For example, highly LREE-enriched REE + Y patterns comparable to the ORB waters were documented in marine stromatolites of a restricted ca. 2.82 Ga epi-continental basin, and interpreted by Kamber et al. (2004) to reflect highly localized terrestrial REE + Y input from tonalitic gneisses.

5.3 Terrestrial Fractionation of Zr/Hf and Nb/Ta

Significant Zr-Hf fractionation and a general trend towards superchondritic Zr/Hf ratios is evident across rivers of variable ionic strength, pH, and Eh (Godfrey et al. 1996; Inguaggiato et al. 2015; Zuddas et al. 2017; Censi et al. 2018; Zuddas et al. 2018). Ionic strength and the composition of colloids appear to play major roles in Zr-Hf cycling within the dissolved load of rivers (Censi et al. 2018). Both Zr and Hf are predicted to occur primarily as M(OH)5 with minor M(OH)04 hydroxyl species (Turner et al. 1981; Byrne 2002) in near-neutral waters, with Hf > Zr adsorption to positively charged (oxy)(hydr)oxide surfaces seemingly a key process generating superchondritic Zr/Hf in both terrestrial and marine environments (Godfrey et al. 1996; Firdaus et al. 2008, 2011; Schmidt et al. 2014; Censi et al. 2018; Zuddas et al. 2018).

The Zr/Hf ratios of SLRS-6 and the Ottawa River T-C transect scatter within the range of upper crustal rocks despite the downstream decrease in Zr-Hf abundances. The elements are transferred from their crustal source into the dissolved load (assumedly affiliated with colloids) and transported in the Ottawa River without significant fractionation. The indistinguishable Zr/Hf ratio of SLRS-6 from the T-C transect also suggests the inherited crustal signature is consistent across different levels of filtration (< 0.45 µm for this study, < 0.2 µm for SLRS-6) and carried to the southernmost ORB. This could imply a limited chemical reactivity with riverine particulates and physical colloid removal (e.g. colloid aggregation), as suggested by Godfrey et al. (2008) for low salinity samples of a Hudson River estuary transect. In contrast to the Ottawa River, the subcrustal Zr/Hf ratios of the 4 ORB tributaries indicate preferential Hf > Zr transfer into the dissolved load via weathering or some colloid-particulate interaction that favours Hf > Zr aqueous stability. This direction of fractionation is similar to a subset of low ionic strength waters studied by Censi et al. (2018) and highly acidic waters studied by Inguaggiato et al. (2015). In contrast to these studies, assessing the role of specific ligands or colloid types was not a focus here. Nevertheless, it is noted that independent of the different Zr/Hf signatures in ORB tributaries (with lower Zr-Hf abundances) they must ultimately be overwhelmed by the signature of the trunk Ottawa River downstream. There are no Zr-Hf data available to trace the coupled pathway of these elements further into the St. Lawrence River and St. Lawrence estuary, but Gobeil et al. (2005) provide evidence that the natural Zr flux from the wider Great Lakes catchment is retained with negligible anthropogenic modification in the St. Lawrence River downstream of Montreal.

A comparison of the ORB Zr-Hf geochemistry to seawater and other freshwater is presented in Fig. 11. In general, river waters have overlapping or higher Zr abundances than deep seawater but more limited Zr/Hf fractionation from crustal values. While the new results from this study do not provide new details on specific mechanisms of terrestrial Zr/Hf fractionation, they emphasize the importance of varying Zr/Hf signatures across rivers with differing chemistry and catchment geology. Specifically, the close overlap of the ORB waters with those of Godfrey et al. (1996) and Godfrey et al. (2008) suggest riverine fluxes from silicate-dominated catchments with low ionic strength may be prone to have minimally fractionated Zr/Hf (relative to crustal rocks) compared to higher ionic strength catchments dominated by carbonates and evaporites (Zuddas et al. 2017; Censi et al. 2018; Zuddas et al. 2018). The study of Tepe and Bau (2014) demonstrated that enhanced fine particulate and colloid input to rivers could play a dominating role in the dissolved (< 0.45 µm) input to local marine waters. It is likely a combination of dissolved riverine fluxes with varying Zr-Hf abundances and Zr/Hf ratios, Zr-Hf released from atmospheric dust particles (Censi et al. 2019), and coastal sediment Zr-Hf cycling that contribute to the significant Zr/Hf variability observed in shallow ocean waters (Firdaus et al. 2011; Niu 2012) and thus further work in each area is still needed.

Fig. 11
figure 11

Compilation of seawater and freshwater dissolved Zr–Hf–Nb–Ta data compared to the ORB waters in a plot of Zr versus Zr/Hf (a) and Nb versus Nb/Ta (b). Seawater is divided into shallow and deep at a depth of 1000 m as per Niu (2012). Compiled seawater data from: NE Atlantic Ocean (Godfrey et al. 1996); WN Pacific Ocean (Firdaus et al. 2008); Pacific Ocean (Firdaus et al. 2018); N Atlantic Ocean (Firdaus et al. 2018); NE Indian Ocean and surrounding seas (Firdaus et al. 2018). Compiled freshwater data for Zr-Hf (Godfrey et al. 1996, Firdaus et al. 2008, Censi et al. 2018) show higher Zr and lower Zr/Hf relative to seawater. Freshwater Nb–Ta data are highly limited with only the Uji River (Firdaus et al. 2008) and SLRS-4 (Filella and Rodushkin 2018) available to compare to the SLRS-6 and ORB data from this study

Insight into the Nb–Ta geochemistry of rivers is restricted to a handful of data (Filella 2017; Filella and Rodushkin 2018). Most Nb/Ta studies on seawater to date draw a constraint for riverine Nb/Ta ratios from Firdaus et al. (2008) based on only SLRS-3 data (Nb: 2.9; Ta: 0.25; Nb/Ta: 11.4) and the Uji River (Nb: 7.0; Ta: 0.69; Nb/Ta: 10.1). However, as illustrated in Sect. 4.6.2, even the SLRS series CRM Nb/Ta ratios are poorly constrained with highly variable inter-study results. As such, the Nb–Ta geochemistry of rivers and specific controls of certain ligands or colloid types across different compositions of water is largely unconstrained compared to Zr-Hf. Both Nb and Ta are predicted to occur primarily as M(OH)05 hydroxyl species in near-neutral waters (Turner et al. 1981; Byrne 2002; Koschinsky and Hein 2003). To explain oceanic Nb/Ta fractionation, Schmidt et al. (2014) suggested that the greater stability of hydrolyzed Nb relative to Ta contributes to generating supercrustal ratios.

The Nb/Ta ratios of SLRS-6 and the Ottawa River T-C transect are consistently supercrustal and increase downstream with decreasing Nb–Ta abundances. Both observations are surprising and contrast with Zr-Hf in illustrating Ta > Nb removal throughout the catchment. This is also surprising in view of available data indicating that Nb/Ta fractionation is more subdued relative to Zr/Hf in the oceans (Firdaus et al. 2011, 2018). A coupling of Nb/Ta with changes in HFSE abundances and Cen/Cen* provides supporting evidence that fractionation is likely tied in some way to colloid removal process in the ORB, but further detail on mechanisms is not speculated on here.

A comparison of ORB Nb–Ta geochemistry to seawater and other freshwater is presented in Fig. 11. Freshwaters are generally higher in Nb abundance than shallow and deep seawater, and earlier observations suggested near-crustal Nb/Ta ratios dominate in rivers, implying that most Nb/Ta fractionation occurs in shallow seawater. However, the new ORB data from this study and the recently published Nb/Ta ratio for SLRS-4 (26.4 ± 5.9) from Filella and Rodushkin (2018) overlap with part of the supercrustal Nb/Ta range observed for seawater. Thus, the results of this study tentatively suggest that the very low budget of dissolved (colloidal) Nb–Ta transferred to the hydrosphere begins to fractionate in some terrestrial systems prior to the oceans. Negligible Nb/Ta fractionation has been reported for the residual Nb–Ta in a saprolitic weathering profile (Babechuk et al. 2015), but, overall, most steps in the source-to-sink aqueous Nb–Ta cycle are very poorly understood. Better constraints on riverine Nb–Ta fractionation and effects associated with coastal colloid-particulate processes and remineralization are needed to derive a more complete understanding of oceanic values and the potential utility of Nb/Ta ratios as a marine biogeochemical tracer (Firdaus et al. 2011, 2018).

5.4 Terrestrial Fractionation of Th/U and Mo/W

Transport of dissolved U in oxygenated and near-neutral aqueous systems is primarily related to the oxidation of U(IV) to U(VI) and its solubilization as the uranyl cation (UO22+) or a uranyl complex (e.g. UO2(CO3)2−2 ) and potentially aided by complexation with organic ligands (Halbach et al. 1980; Lenhart et al. 2000) and siderophores (Kraemer et al. 2015). In contrast to U, Th is redox-insensitive and significantly less soluble, ultimately concentrating in the particulate/colloidal fraction of oxygenated waters (Moulin and Ouzounian 1992; Porcelli et al. 2001; Pokrovsky et al. 2010). These differences result in Th-U fractionation throughout soils and rivers (Mathieu et al. 1995; Chabaux et al. 2003; Suhr et al. 2018; Yu et al. 2019) relative to upper crustal rocks with a Th/U ratio ~ 4 (Kamber et al. 2005; Rudnick and Gao 2014). This fractionation leads to generally subcrustal Th/U ratios in river water, as evident from the global riverine Th/U data compiled by Gaillardet et al. (2014) ranging from 0.01 to 6.23 with a mean of 1.38 ± 3.26 (n = 24).

The consistently subcrustal Th/U ratios in the dissolved load of the ORB waters indicate the preferential chemical weathering transfer and solubilization of U > Th expected for oxygenated conditions and widely observed in other catchments (e.g. Yu et al. 2019). The selectively higher Th/U ratios of the two tributaries draining from the east into the Ottawa River (Rivière Noire and Coulonge) suggest a lithological control releasing less U or a catchment processes leading to enhanced U retention (e.g. more reductive scavenging in organic-rich waters/sediment) in comparison to other ORB waters. The similar Th/U ratio of SLRS-6 to the new OR samples (Fig. 9) suggest that signatures inherited from silicate-dominated parts of the catchment are not significantly influenced downstream by U-bearing Paleozoic sedimentary rocks or agricultural contaminants that could produce even lower subcrustal Th/U ratios.

Both Mo and W are redox-sensitive with a solubility aided by oxyanion (MoO42−, WO42−) formation during weathering and aqueous transport in oxygenated and near-neutral environments (Turner et al. 1981). The Mo/W of the upper continental crust composite of Rudnick and Gao (2014) is 0.58 and alluvial sediment of the Murray-Darling Basin, Australia (Marx and Kamber 2010) is slightly lower with a Mo/W mean of 0.28 ± 0.23 (n = 105). The Mo/W ratio of seawater is significantly higher at ~ 950–1150 and rather constant with depth in oxygenated waters due to their near-conservative marine behaviour and long residence time (Sohrin et al. 1987; Firdaus et al. 2008). The much higher Mo/W of seawater than upper crust is attributed primarily to the preferential scavenging of W > Mo on Fe–Mn (oxy)(hydr)oxide surfaces, and especially those of Fe (oxy)(hydr)oxides, under oxic-hypoxic conditions (Takematsu et al. 1990; Hein et al. 2003; Kashiwabara et al. 2013, 2017; Dellwig et al. 2019). However, Mo > W removal to sediment also occurs in sulfidic environments (Bauer et al. 2017). Prior to the oceans, decoupling of Mo-W can be traced across estuaries where Mo behaves more conservatively than W (van der Sloot et al. 1985; Bauer et al. 2018) as well as into terrestrial environments as indicated by river waters with Mo/W ratios ranging from 0.3 to 68 (Mohajerin et al. 2016; Bauer et al. 2018). Terrestrial fractionation appears to begin in weathering profiles, but relative Mo-W mobility is probably controlled initially by protolith mineralogy (Greaney et al. 2018) as evident by preferential W > Mo release recorded in a basaltic saprolite (Babechuk et al. 2015). The details of pedogenic Mo-W fractionation are poorly understood.

The consistently supercrustal Mo/W ratios of the ORB waters (12.0–74.5) cover a similar range to other river waters measured to date (Mohajerin et al. 2016; Bauer et al. 2018). The entire Ottawa River T-C transect has supercrustal Mo/W ratios with minimal downstream variation, suggesting that the signature is inherited primarily from processes in aquatic environments nearer the source, or in soils. Assuming the bedrock in the ORB has an average Mo/W ratio of < 1 similar to upper crustal estimates, the ORB water ratios indicate either incongruent weathering favouring Mo release or W > Mo adsorption effects on Fe (oxy)(hydr)oxide surfaces (comparable to marine environments). However, more organic-rich or lower pH conditions could also favour increased Mo/W ratios (Bednar et al. 2009). Some combination of these processes appears to produce the more extreme Mo/W ratios locally in the tributaries. Without further work, the mechanisms of Mo-W fractionation are speculative. Nevertheless, the results suggest that Mo/W ratios could emerge as a reasonable tracer of terrestrial (oxy)(hydr)roxide development in a similar manner to Ce anomalies. The minor downstream variations of Mo/W in the Ottawa River could be related to changing tributary inputs with variable Mo/W ratios; the inflection point where Mo/W ratios start to increase downstream is beyond the area where the Rivière Noire and Rivière Coulonge with higher Mo/W meet the Ottawa River. However, riverine colloid-particulate effects cannot be ruled out in generating some of the Mo/W variation in ORB waters, since W has been documented to have a high affinity for Fe-colloids in rivers and groundwater in contrast to Mo being present as a truly dissolved species (Pokrovsky and Schott 2002; Pokrovsky et al. 2006; Johannesson et al. 2013).

Independent of W, Mo is of notable interest in the ORB due to the unusual Mo stable isotope composition reported for the Ottawa River via the SLRS CRM (Archer and Vance 2008). On an array of river waters from across the globe, the Ottawa River forms an end-member with low Mo abundance and significant enrichment in heavy Mo isotopes. The origin of this signature was interpreted by Archer and Vance (2008) to indicate the preferential retention of light Mo isotopes in soils (King et al. 2018). However, a detailed assessment of other potential effects such as inheritance from atmospheric particulates (King et al. 2016) or anthropogenic/natural input from ore deposits in the northern ORB has yet to be undertaken.

6 Conclusions

This study describes a simple and rapid ICP-MS workflow capable of producing high-precision natural water data for ≥ 38 analytes (Ca, Mg, Na, K, and 34 trace elements) down to pg g−1 levels from sample volumes ≤ 9 mL. This direct analysis method minimizes preparation time, costs, and sample handling and associated blank compared to pre-concentration strategies (Bau and Dulski 1996; Bayon et al. 2011; Fisher and Kara 2016; Hoang et al. 2019), but it requires access to clean laboratory space to reduce blanks in reagents and containers used in sampling and analysis. The final analyte list can be extended, or modified in terms of acquisition parameters, depending on the volume of sample and instrumentation available. The method is customizable to even lower sample volumes than the lowest (~ 2 mL) used here (at a compromise between uptake rate, number of analytes, and acquisition time) and adaptable to other types of ICP-MS and sample introduction strategies offering higher mass resolution or enhanced interference removal (e.g. sector field ICP-MS or ICP-MS/MS instruments). The method is also appropriate for saline samples following pre-dilution to minimize matrix ion effects during analysis (see Lawrence and Kamber 2006, 2007). All of these factors make the method attractive for wide regional sampling programs (e.g. Leybourne and Johannesson 2008), low-volume sample applications (e.g. extracted pore waters, soil moisture, plant fluids, etc.), or studies requiring the bulk of the sample to be preserved for other analyses (e.g. stable metal isotope ratio analysis). However, detection limits for the lowest abundance elements (e.g. the HFSE Ta, Hf, W) are liable to be reached for very trace element-poor natural waters or when applying pre-dilution.

A compilation of NRC-CNRC SLRS-6 river water CRM trace element data (63 measurements from 7 independent bottles) measured with the ICP-MS method across two different facilities provides constraints on method precision and accuracy, as well as new recommended SLRS-6 information values. Inter-bottle variations were negligible within analytical uncertainty. Precision was demonstrated to be excellent for the vast majority of analytes (e.g. ≤ 5% 2RSD for Na, Mg, K, Ca, Li, Rb, Sr, Zr, Nb, Mo, Cs, Ba, REE + Y, Hf, W, U). Of the elements measured, 11 are certified by NRC-CNRC, and agreement was excellent for 9 elements (bias < 5%: Na, Mg, K, Ca, Sr, Sb, Ba, Pb, U), good for 1 element (bias 5–10%: Mo), and poor for 1 element (bias > 10%: Cd). For 27 uncertified elements, comparison with literature values from Yeghicheyan et al. (2019) yields excellent to good agreement (< 10% bias) with the exception of some of the lowest abundance elements (Tl, Th, Hf, W). For several of the HFSE, the results of this study are generally lower than in Yeghicheyan et al. (2019) and agree much closer with reported abundances from earlier generations of the SLRS CRM. The lack of certified HFSE abundances (notably Nb–Ta-Zr-Hf-W) in the SLRS series and other natural water CRM remains a major barrier in method development and inter-laboratory comparison. Further, the long-term aqueous stability of some HFSE (e.g. Ta, Th) in distributed CRM bottles is unclear. To address both of these issues, it is recommended that an analytical round-robin involving facilities specialized in low-level HFSE measurements be organized shortly after a new CRM is made available, and for each facility to subsequently run long-term tests of HFSE stability.

The first ultra-trace element characterization of the dissolved (< 0.45 µm) load of ORB waters (14 samples from the Ottawa River, 11 from tributaries, and 4 from lakes/ponds) was undertaken, with focus on areas at and upstream of the SLRS CRM collection sites. These data provide context to the signatures measured for decades in the SLRS CRM series. Beyond this, the data offer a snapshot of riverine signatures generated in a catchment draining predominantly Precambrian felsic igneous/metamorphic rocks. The catchment geology produces an atypically LREE > HREE-enriched upper continental crust-normalized REE + Y pattern in the dissolved load of the ORB river waters. A downstream decrease in REE + Y abundance points to either mixing/dilution effects or colloid removal in the Ottawa River that ultimately results in more negative Eu and Ce anomalies but otherwise does not significantly change the REE + Y pattern. The HFSE (Nb–Ta-Zr-Hf–Th-U–Mo-W) show variable fractionation relative to estimated bulk compositions of the upper continental crust, in accordance with their relative solubility, riverine particle reactivity, and presumed affiliation with colloidal material. Zirconium and Hf abundances decrease downstream in the Ottawa River, yet crustal Zr/Hf ratios are retained and indicate weathering and colloidal transport/removal processes do not result in significant Zr/Hf fractionation, similar to some previous findings (Godfrey et al. 1996; Firdaus et al. 2008; Godfrey et al. 2008) but contrasting with others from higher ionic strength waters (Zuddas et al. 2017; Censi et al. 2018; Zuddas et al. 2018). Niobium and Ta abundances decrease downstream in the Ottawa River, but with consistently supercrustal and increasing Nb/Ta ratios. These results appear to capture hitherto undescribed terrestrial Nb/Ta fractionation. However, with a highly limited availability of Nb/Ta ratios from other rivers for comparison, the full extent and environmental importance of terrestrial Nb/Ta fractionation remains open (Filella 2017). The Th/U and Mo/W ratios of ORB waters are consistent with preferential transfer and solubilization of U > Th and Mo > W. These data provide further evidence for terrestrial Mo/W fractionation, which is still poorly understood, and suggests it could indirectly track soil (oxy)(hydr)oxide development during oxidative chemical weathering akin to Ce anomalies, or some localized aquatic particulate W > Mo removal upstream of sampling.

The combination of multiple REE + Y and HFSE fingerprints reveal different aspects of the ORB catchment weathering and hydrology. These results illustrate one of the applied advantages of the ICP-MS method – the ability to measure a full suite of elements at ultra-trace (< ng g−1) levels in a single, rapid analysis. Without the assessment of particulate/colloid fraction geochemistry or other hydrological parameters (e.g. anion budgets, mineral saturation modelling) for this study, the results are designed primarily to layout the baseline ultra-trace element geochemistry of the ORB for more in-depth and targeted studies in the future. Future studies are recommended to better define processes in ORB soils, include a more complete characterization of the upper ORB (notably east of Lake Timiskaming and closer to the Ottawa River head waters), and revisit ultra-trace element signatures with apportioning between colloidal, suspended particulate, and truly dissolved loads.