Detecting trends in freshwater trace element concentrations: methodological issues and data treatment

Anthropogenic activity has significantly perturbed the Earth's natural biogeochemical cycles, resulting in a significant translocation of many chemical elements from the Earth's crust to human-built objects and to emissions to the environmental compartments. In some environmental media, changes in trace element (TE) concentration have been related to their use. This is the case for sediments, mosses, ombrotrophic peatbogs, or ice cores. It is often assumed that increased human use of chemical elements should also lead to increasing concentrations in surface continental waters (far from direct pollution sources), but this is not clear-cut. Temporal trends in concentrations of major ions have been studied, mostly in the context of acid rain recovery (Garmo et al. 2014; Vuorenmaa et al. 2017) or of organic carbon evolution (Filella & Rodríguez-Murillo 2014), but studies on temporal trends of TE in surface waters are comparatively very scarce (Huser et al. 2011; Todd et al. 2012). The main reason is probably that there are few good, analytically traceable, long-term data series of TE concentrations – like those required for long-term trends to emerge over temporal variations due to occasional changes, seasonality and analytical errors. In particular, TEs have suffered from a lack of appropriate analytical methods until relatively recently (e.g., quadrupole-based inductively coupled plasma–mass spectrometry (ICP-MS) was introduced commercially in 1983) and has benefited from improved trace metal clean procedures introduced in the early 1990s.

Although often just visual (Hatje et al. 2016) or linear correlation methods are applied, rigorous work requires the use of non-parametric methods and the consideration of the possible existence of serial dependence and long-term persistence (LTP). The purpose of the present study is to explore methodological and data treatment issues when attempting to rigorously determine temporal trends in freshwater TE concentrations. For this purpose, we have used a temporally extended, traceable dataset of waters from Lake Geneva, Switzerland. TE concentrations have been collected since 1994 in this lake as part of a water-quality monitoring program set up by the Service de l'Ecologie de l'Eau (SECOE). Lake Geneva is a particularly interesting system because it is a high volume water body with an elevated water residence time, situated in an anthropised environment but not subjected to large direct pollution sources.

Seven elements were chosen for this study (Cd, Co, Gd, Mo, Pb, Sb, Sr). They cover different types of chemical behaviour (Supplementary Information, Figure SI1) and they also differ with respect to the type and intensity of their use (Supplementary Information, Figure SI2). In particular, Gd and Pb are interesting reference points. The main use of Pb (i.e., leaded gasoline) was completely phased out in all countries over a relatively short period of time that is well documented (Filella 2017). A noticeable, widespread increase of Gd concentration in surface waters in the past decade has been world-wide observed, likely due to the element's use in magnetic resonance imaging contrast agents (Kulaksiz & Bau 2011; Hatje et al. 2016).

Lake Geneva lies at an altitude of 372 m and has a surface area of 580 km² (Supplementary Information, Figure SI3). It is the biggest body of freshwater in Western Europe and supplies drinking water for more than 800,000 people. Its average water residence time is 11.3 years. Further details can be found in Table SI1 in the Supplementary Information file.

Monthly measurement of TE concentrations in point G3 (located at 6.2197° E/46.2994° N; water depth: 72 m) started in December 1995. Measurements are taken at ten depths from 0 to 70 m. Water samples are collected with a PVC Kemmerer sampling bottle lowered in the water from a research vessel. Samples for metal analysis are filtered through 0.45 μm pore-size filters (Whatman^® Nuclepore), acidified to pH 1 with Suprapur nitric acid and kept at 2 °C(±2) until analysis. Thus, this study is based on so-called ‘dissolved’ concentrations. For some elements, that are present mostly dissolved in freshwaters, such as antimony (Filella et al. 2009), this is not relevant, but for others, such as lead, that have been shown to be present in variable but never negligible amounts in the particulate fraction (Town & Filella 2002), this fact needs to be taken into consideration when interpreting the results.

Concentrations of TE were measured by ICP-MS on a Fisons PQ2+ before 2006 and on a ThermoFisher X7 II from 2006 to 2015, both in their standard configurations at the Service de l'Ecologie de l'Eau (SECOE) of Canton Geneva, Switzerland. The analytical accuracy was followed by analysing the certified reference materials (CRM) SLRS-4 and SLRS-5 (National Research Council Canada, https://www.nrc-cnrc.gc.ca/eng/solutions/advisory/crm/list_product.html). Detection limits (μg L⁻¹) are: Cd, 0.001; Co, 0.01; Gd, 0.0005; Mo, 0.05; Pb, 0.001; Sb, 0.002; Sr, 0.2. The analytical laboratory is accredited for the analysis of the trace elements in freshwaters (accreditation number: STS245; ISO norm 17025).

Data are shown in Figure SI4 in the Supplementary Information file. Temporal trends have been evaluated for all data (all depths together; this dataset will be named All data), integrated element load (IEL) and TE concentrations at each depth. IEL is calculated by: (i) dividing the water column into layers, each layer containing all the concentration data measured at equally spaced points; (ii) averaging TE concentrations in each layer to obtain the mass of TE in a column of 1 dm² of section area and the height of the corresponding layer; and (iii) adding all TE masses from surface to bottom. Integrated TE is then defined as the above mass divided by the lake depth in dm; this represents the weighted TE concentration at the sampling point.

Non-parametric methods are required to analyse temporal trends of TE concentrations in water because data are usually non-normally distributed (Reimann & Filzmoser 2000). Thus, temporal trends have been studied by applying Mann–Kendall (MK) or seasonal Mann–Kendall (SMK) methods (Hirsch & Slack 1984; Helsel & Hirsch 2002), depending on the presence or not of seasonality in the data. Seasonality is expected in lakes, particularly for biogenic elements (Sr, Mo, Co). Here, it has been assessed using the Kruskall–Wallis test of equal medians for the data classified by depth and month. The test was applied to data at 0 and 70 m. Only Sr concentrations at 0 m data showed seasonality (p = 1.5 × 10⁻⁶ < 0.05) (p = 0.074 at 70 m depth). Thus, MK was applied for all TE except Sr. IEL data were also assessed for seasonality, which again was found only for Sr (p = 0.019).

The statistical significance of the trends is obtained in MK methods by applying the Z-statistic test of the sum of signs of the differences between every pair of values; Z shows a normal distribution and therefore the statistical significance of the temporal trend can be evaluated by the p value. However, blind use of p values is being increasingly challenged, particularly in the case of tiny effects (Siontis & Ioannidis 2011). We therefore show all trends, whether significant or not, always accompanied by their corresponding Z value, which is directly related to p.

The average change over time of TE concentrations was calculated using Sen's slopes (Sen 1968). A Sen's slope is the median of the slopes between all possible pairs of temporal data (N(N−1)/2), N being the number of pieces of data in the series.

Measured values are very close to detection limits for some elements (Cd, Gd and Pb) and their data series include some censored data (Supplementary Information, Table SI2). Data have been considered censored when equal to or less than the minimum positive value in each data series. When possible, Sen's slopes have been calculated using the non-parametric method of Akritas et al. (1995) for data at each depth to take the effect of censoring into account.

Median TE concentration profile values have been calculated for the period 2006–2015. Significant differences between Sen's slopes at different depths were assessed by the bootstrap technique following Wilcox (2012). In brief, 10,000 samples (i.e., concentration data series) are obtained by resampling – random picking of elements with replacement – from each time series at each depth, giving 1,000 bootstrap Sen's slopes (SS_b). Sen's slope at another depth j (SS_j) is then compared with Sen's slope of a depth i (SS_i) by determining whether SS_j is within 95% of SS_b.i. If not, the null hypothesis H₀ = ‘equal slopes’ can be rejected at the 1−α confidence level (i.e., SS_j and SS_i are different). The same technique was applied to estimate the significance of differences between median concentrations at the ten depths.

Hamed and Rao's method (Hamed & Rao 1998) was used to correct MK when there was serial correlation and the application of the method was possible.

Programs for MK and SMK tests are from www.mathworks.com/matlabcentral/fx_files/22389/7/sktt.m, www.mathworks.com/matlabcentral/fileexchange/11190-mann-kendalltau-b-with-sens-method-enhanced/content/ktaub.m adapted by us. Kruskall–Wallis is from PAST (Hammer et al. 2001). Bootstrap calculations have been performed with built-in Matlab functions. Hamed and Rao's method was implemented in Matlab using the ‘Mann-Kendall_Modified_Fix.m’ function, developed by Manjurul Hussain Shourov (https://www.researchgate.net/post/Does_anyone_know_the_package_for_Modified_mann_kendall_trend_test_in_R_software) that corrects a bug in http://www.mathworks.com/matlabcentral/fileexchange/25533-mann-kendall-modified-test. Sen's slope calculations and bootstrap for series with censored data have been calculated using R-3.4.1. For slopes, the ‘cenken’ function in the ‘NADA’ package was used (Helsel 2012).

TE concentrations for Cd, Co, Mo, Pb, Sb, Sr (1995–2015) and Gd (2006–2015) at all depths in Lake Geneva (G3) (All data) are shown in Figure 1. The dramatic improvement in the precision of the measurements due to the change of the ICP-MS instrument in 2006 is readily visible. This effect is noticeable for all TEs considered and is particularly strong for Cd and Sb.

Before any data were treated, the autocorrelation structure of temporal series of IEL data was evaluated. Autocorrelation functions (Supplementary Information, Figure SI5) show that, except for Cd, significant serial correlation exists for all elements and, consequently, the Hamed and Rao method to evaluate the significance of the temporal trends needs to be applied.

No visually apparent change of tendency is seen along the whole studied period in All data concentrations of three elements (Mo, Sb and Sr) (Figure 2). The change of ICP-MS apparatus in 2006 clearly introduced an offset in the data values for Mo and Sb. Data for these elements were treated separately before and after this year to check for possible spurious effects of the instrument change on the trends observed (Table 1). Although no shift is visible for Sr, data for this element were also treated separately. No significant trends exist for Sr. Increases were observed for All data concentrations and IEL for Mo and Sb in all periods, but they became non-significant for IEL before 2006 when Hamed and Rao's corrected Z values are considered. Gd concentration increases significantly (2006–2015), both considering All data and IEL.

Changes in direction of the trend during the sampling period are apparently visible for Cd, Co and Pb (Figure 1). Their concentrations appear to decrease in the first half of the period considered and are stationary or increase in the second half. Since the MK method requires a monotonic trend (otherwise the trend could become dependent on the initial and final dates of the time series (Stevenson et al. 2012)), the two data periods were treated separately. Turning years were 2001 for Pb and 2006 for Cd and Co. All data concentrations and IEL of Cd and Co do not change significantly before the turning year and increase significantly afterwards. However, the IEL increasing trend for Co becomes non-significant when autocorrelation is considered. All data concentrations and IEL for Pb decrease in all periods, although the change in IEL is not significant once the Hamed and Rao correction is applied.

Net and percentage changes in IELs can be found in Table 1 for the period 2006–2015.

Lake data at different depths were evaluated only after 2006 in order to avoid noisy data. Median concentration profiles of TEs in the lake (2006–2015) are shown in Figure 2. Results of the bootstrap data treatment are shown in Table SI3 in the Supplementary Information file. Co, Gd and Mo show no statistically significant changes of their median concentrations with depth while Pb concentrations both at the surface and at the bottom are significantly higher than at other depths. Sb also shows differences (higher concentration) at 0 m compared with 20 and 70 m. Sr median concentration values are higher below 20 m. Cd concentrations are lower at 2.5 m but significantly only when compared with a few other depths (20, 50 and 70 m).

Sen's slopes at all depths (2006–2015) are also shown in Figure 2. Numerical values are given in Table SI4 in the Supplementary Information file. All are positive, except for Pb at the surface, and statistically significant at all depths in the case of Cd (except at one depth), Gd, Mo and Sb. For Sr, they are only significant at certain depths. Differences in Sen's slope values between depths, evaluated by bootstrap, are shown in Table SI5 in the Supplementary Information file. Co and Sr show no variation at different depths. Mo, Pb and Sb have lower values at the surface compared to all depths in the case of Pb and to some deeper levels in the case of Mo (7.5, 15, 20 and 30 m) and Sb (15, 20, 50 and 70 m). Some elements show different Sen's slopes at the bottom, smaller in the case of Mo (compared to slopes at 7.5 and 20 m) and bigger for Sb (compared to slopes at 0, 5 and 30 m).

Table 2 summarises the results obtained for the various elements, periods and parameters considered. It reflects the complexity of the data treatment, particularly the difficulty of assigning statistical significance to weak trends. The work carried out has made it possible to establish the temporal trends for TEs in Lake Geneva, but also to gain a better understanding of the methodological caveats that apply to temporal trend evaluation of TEs in freshwaters. These are discussed below.

Data in trend studies have usually been obtained from monitoring programs rather than from research-oriented projects, which implies different objectives and data quality standards. Nevertheless, our study shows that temporal trends can be followed using monitoring data, provided that it is possible to have direct access to the raw data and protocols used (i.e., changes of instrument, calibration, use of CRMs, data pre-treatment). However, detailed information on quality assurance (QA) all through the years is not available even in a ‘favourable’ case like ours. Unfortunately, this is a problem inherent to any data trend study based on monitoring data.

How long does a data series need to be for a temporal trend to be detected? This depends on two factors: the strength of the change and the quality of the data. It can be roughly evaluated by calculating the time of emergence (ToE) parameter, which gives the point in time when a signal finally emerges from the background noise of natural variability (Keller et al. 2014). Values for ToE of the elements studied (Table 3), calculated before and after 2006, clearly show that where data are noisier (i.e., before 2006), longer time periods are needed for the trend to be detectable. In the case of Pb and Sr, ToEs are higher in the second period because trends are stronger in the first period and compensate for the greater noise. As expected, elements with significant changes have shorter ToEs, closer to ten years.

Pb, Cd and Gd data have a sizeable fraction of censored data (Table SI2). Simply substituting a fabricated value (zero, detection limits, or another value) for censored data can lead to wrong results in statistical data (mean, standard deviation), hypothesis tests and regression models (Helsel 2012). Median concentrations are not affected by left-censored data (i.e., those less than a limit) if censoring is less than 50% of the data. Although trends are most effectively detected in uncensored data as compared to censored data, the existence of censored data will not give any spurious significant trend (i.e., its effect is eliminating information). This means that statistical significance in MK is not affected by censoring. However, it has an effect on the value of Sen's slopes (Table SI2). Those slopes have thus been calculated using the non-parametric method of Akritas et al. (1995) for data at each depth, but calculations could not be performed for All data because the method does not work when tied values are present (in our case, there are ten concentration values, i.e., ties at each time).

Another effect of data quality is related to the precision of the measurements that, when close to the detection limit, results in very few significant figures in the readings. Here lies the reason behind the zero Z values for Cd and Co data before 2006 (remember that MK is a ranking method) and the ‘ugly’ aspect of Gd in Figure 1. A final aspect that has, to our knowledge, never been explored in concentration trend studies, is the effect on the measured concentration values of potential drifts of the instruments over the years. We have tried to follow potential drifts by applying the same data treatment to the CRM concentrations measured (Table SI6), with no conclusive results (see discussion in the Supplementary Information).

Finding statistically significant, deterministic temporal trends in time series is complicated by two features of time series: serial correlation and LTP. Hydrological time series, such as those we are dealing with, very often show positive serial correlation. Serially correlated data without any deterministic trend tend to stick above or below the average value of the series more frequently than would be the case for a purely random process. This increases the likelihood of identifying a significant trend where there is none. Hamed and Rao's method (Hamed & Rao 1998) was used to correct MK when there was serial correlation, but it could only be applied to IEL data and not to data from all depths lumped together because the method does not allow for ties (more than one value at each time). As expected (Cox & Stuart 1955), the application of the corresponding correction unequivocally leads to a decrease in the significance of the observed trends when it could be applied (i.e., IEL).

A ‘persistent’ time series shows long periods of values above or below the mean, even if the mean does not change over the long term (Beran 1994). Unfortunately, deciding whether LTP is present and accounting for it requires data series longer than the ones we are studying here (O'Connell et al. 2016). We therefore did not correct our results for LTP. A deeper discussion on these problems can be found in the Supplementary Information file (Section SI1).

Trend studies commonly attempt to ascribe observed effects to well-defined, physical causes. This is usually done by establishing correlations, and backing them up by informally including other lines of evidence from the literature in the Discussion section (Downes et al. 2002). Apart from the well-known fact that correlation is not causation, in practice, this approach is not straightforward. Only in rare cases is the trend strong enough and due to a few, well-defined, drivers (e.g., recovery from acid rain); generally, observed effects stem from the interaction of a variety of interrelated causes that are not always acting in the same direction, often leading to ‘weak’ causation effects. This type of exercise or the application of more sophisticated approaches (Granger 1980; Kleinberg 2012) is entirely outside of the scope of this paper. It is interesting however to discuss what effects can be reasonably expected and whether the results obtained fit with existing knowledge.

Temporal changes in lake TE concentrations will depend on temporal changes in lake TE inputs (i.e., direct input from waste water treatment plants, direct atmospheric deposition on the lake, input from the lake catchment) and/or temporal changes of internal physical and chemical processes taking place in the water column and sediments. Obviously, not every change in an element input will result in a statistically significant increase in concentrations because, on the one hand, increasing input load might be compensated for by internal lake processes while, on the other, the input itself must be (i) big enough compared to other sources, (ii) big enough compared to lake concentrations and (iii) present in a ‘convenient’ form. Examples of these conditions are found in our results.

Big enough compared to other sources: There are 171 WWTPs in the Lake Geneva watershed, which is home to a 3,009,830 population equivalent (Table SI2). A recent quantitative assessment of TE flows from WWTP in Switzerland (Vriens et al. 2017) showed that the contribution of WWTP effluents to river flow is very low for most of the elements we are concerned with here: Cd (3%), Co (8%), Mo (5%), Pb (1%) and Sr (2%). The WWTP contribution is a bit higher for Sb (16%) and huge for Gd (83%). In the case of Gd, it is widely accepted that this is due to its medical use (Kuroda et al. 2016) and, for Sb, it is probably the result of its widespread use and dissipation in the final products (Turner & Filella 2017). Therefore, the only case where the observed temporal trend can be unequivocally attributed to the input from WWTP is Gd.

Big enough compared to lake concentrations: Direct deposition of TE on the lake will have a measurable effect on surface concentrations only when deposition flux and accumulation thereof are at least an order of magnitude higher than element water concentrations. For instance, if we consider the elements studied here, Sr concentrations in the lake are more than two to five orders of magnitude higher than for the other elements studied, while atmospheric deposition is negligible. It is therefore only to be expected that no effect is observed. Atmospheric deposition would need to be very high indeed for an effect to be seen for this element. In fact, the only element that shows statistically significant higher concentrations at the surface compared to other depths (2006–2015) is Pb, directly pointing to the importance of atmospheric deposition for this element. Atmospheric deposition cannot be excluded for the other elements, but it is not equally translated into higher surface concentrations, except for Sb, albeit less clearly than for Pb. Atmospheric input is probably not ‘seen’ for Mo because of its higher concentration in lake water compared to Pb and Sb, although Mo concentrations in ice cores in Mount Rosa (Swiss–Italian Alps) – a good measure of metal deposition – were similar to Sb ones (at least before 1995) (Barbante et al. 2004).

Present in a ‘convenient’ form: According to these relatively old data in Mount Rosa, Co deposition was about five times Sb deposition and, in spite of both elements having similar concentrations, no change is observed in Co. This could be explained by the higher variability of Co data (Figure 1) or by Co being present mostly in not readily soluble dust, mainly of Saharian origin (Heimbürger et al. 2010). When TE deposition is in aerosol or particulate form, any effect on lake dissolved concentrations will require that the element is readily, or at least partially, dissolved in lake water. This might be limiting not only in the case of Co, but also for Pb in re-suspension of fine lead-contaminated soil dust (Laidlaw et al. 2012) or of the Sb released from brake pads (Iijima et al. 2008).

Temporal changes of TE concentrations in the surface and upper lake layers, as expressed in Sen's slopes, provide additional information on changes in atmospheric deposition. The most extreme effect on surface Sen's slopes is observed in Pb, where the slope changes to a negative, although non-significant, value. Sb and Mo also show (significant) smaller surface slopes. The simplest explanation of smaller Sen's slopes in and near the surface is the decrease of atmospheric deposition on the lake during the period considered. Detailed changes of TE atmospheric deposition on Lake Geneva during 1995–2015 are not known but existing data for Switzerland point to a general decrease. For instance, drastic decreases of Pb (86%) and Cd (66%) deposition and slight decreases for Co have been measured in Swiss mosses (1990–2010) (Harmens et al. 2013). A 23% average decrease was observed for Sb in mosses from seven European countries from 2005 to 2010 (Harmens et al. 2013).

Even if current decreasing atmospheric TE deposition following emission control is widely accepted (Harmens et al. 2010), except for Pb, our results point to an increase, not a decrease in TE lake concentrations. Thus, a further source seems to be at play. A plausible candidate is a change in the release from the lake catchment.

Not all elements are equally released from soils to waters. An indication of their soil mobility is given by the ratio between their concentrations in soils and in waters. These ratios, estimated from median concentration values in European soils (FOREGS) (Salminen et al. 2005) (soil concentrations in mg g⁻¹ and water concentrations in mg L⁻¹), give: Gd (385) > Pb (243) > Co (49) > Cd (15) > Sb (8.6) > Mo (2.8) > Sr (0.81). Very low mobility is thus to be expected for ‘natural’ Gd and Pb as compared to the other elements. Release of TEs deposited from atmospheric deposition or directly added to the soils (e.g., the case of Cd introduced by the use of fertilisers (Roberts 2014)) might differ depending on the form of the deposited element and its interaction with soils but similar behaviour to ‘natural’ TEs is likely. For instance, it is well known that Pb in soils coming from its use in leaded gasoline is extremely insoluble and immobile, with an expected ‘half-life’ of 700 years (Semlali et al. 2004). In any case, changes in TE release from soils, originating either from factors that affect rock weathering (effect of climate change) (Taylor et al. 2012) or from release of TE deposited from atmospheric deposition, are very slow and a delay of years or decades can be expected (Lawlor & Tipping 2003; Huser et al. 2011). Therefore, the origin of the increasing TE concentrations in water (Cd, Mo and Sb) lies very probably in the release of TE previously deposited in soils. A similar effect has been detected in streams and rivers of southern Sweden for Pb (Huser et al. 2011).

Rigorous evaluation of temporal trends of TE in surface waters is not trivial because good, long, traceable concentration data series are required to evaluate whether trends exist. Unfortunately, these data series are scarce, are relatively short because the necessary analytical techniques are quite new, and the quality of the data has until recently often been poor. Moreover, these data have usually been obtained under monitoring programs, not from research-oriented projects, which implies different data quality standards. Nevertheless, studying temporal trends is feasible, provided direct access to raw data is possible and rigorous data treatment methods are applied.

Gd constitutes a good example of how direct input into environmental compartments may have a direct impact, which is easy to detect and explain. However, this type of situation is more the exception than the rule. Neither the often assumed general increase of TE concentrations in environmental compartments linked to the current massive use of the elements, nor the expected decrease, as a result of general, diminishing atmospheric deposition of TE in Europe, are confirmed by this study. This is a nice illustration of the complexity of the response of natural systems to disruption and a warning against over-quick generalisations in environmental sciences.