Cleaner air reveals growing influence of climate on dissolved organic carbon trends in northern headwaters

Water chemistry data were compiled for the period 1990–2016 from regional and national monitoring programs in: Canada (114 sites); Czech Republic (7); Finland (26); the Netherlands (3); Norway (83); Poland (4); Slovakia (12); Sweden (83); the United Kingdom (21); and the United States (75). Records represent mostly headwater lakes and lower-order streams from natural or semi-natural catchments, i.e. largely free of local disturbance, other than forestry and rough grazing practices in sites in the UK [35]. Because we use surface water concentrations of SO₄ and Cl as surrogates for S and Cl deposition, we excluded sites where other processes (e.g. road salt, mine drainage, glacial melting) produce anomalies in those ions. All monitoring programs follow well-established methods and analytical procedures [31, 36]. Sampling frequency ranged from annual to weekly, usually with a higher sampling frequency for streams than for lakes. Previously, temporal patterns in DOC in streams and lakes were found to be similar [37].

For each monitoring site, we calculated annual median concentrations of DOC, as well as several potential chemical drivers of DOC change (SO₄, Cl, pH, nitrate [NO₃], and divalent base cations [∑Ca + Mg]). We used changes in measured surface water concentrations of SO₄ and Cl to represent change in deposition of S and Cl respectively (from both deposition and seasalts). This was because site-specific deposition estimates were generally unavailable for the monitoring sites, modeled estimates often provide wet deposition estimates only, with considerable uncertainty regarding absolute levels, and no estimates were available for the contribution from seasalts. Non-marine concentrations of SO₄ and base cations were also calculated [38] and used as potential chemical explanatory variables.

Gridded daily mean air temperature and precipitation values (0.5° resolution) were downloaded from the Climate Research Unit, University of East Anglia (https://crudata.uea.acuk/cru/data/hrg/; dataset CRU-TS-v3.24.01). Each site was matched to a grid cell (several monitoring sites may be associated with a single grid cell), and daily temperature and precipitation were converted to monthly and annual mean temperature and sum of precipitation. Summer mean temperature and precipitation sum were calculated for July, August and September (JAS).

We applied the Mann–Kendall test [39] to median chemical and mean climate variables to produce estimates of annual trends for DOC and all possible drivers. The magnitudes of the trends were estimated with the Theil–Sen estimator [40]. Trends were calculated for the periods 1990–2016, 1990–2004 and 2002–2016. The periods were determined on the rationale that (a) robust trend analyses require at least ten years of records, but preferably more, (b) a slight overlap prevents that the transition period plays a disproportionate role in the analysis, (c) 1990–2004 was used in Monteith et al [1], thus allowing for a direct comparison with their previously published observations. Sites were included only where trends could be calculated for all three periods of interest, i.e. sites that had (a) data for at least 65% of the total number of years in the period; (b) at least two values in the first four years of the period of interest, (c) at least two values in the last four years of the period of interest. To transform absolute changes in DOC into relative changes in % yr⁻¹, we multiplied Theil–Sen slopes by 100 and divided by the site median concentration for the given period.

Regression analyses were performed using the slope of DOC change (ΔDOC) as the dependent variable, and a range of chemical and climatic variables as candidate predictors, along with interaction variables combining change and state variables (table 1).

Best subsets modeling was performed using the bestglm package in R [41]. Rather than relying on a stepwise approach to adding candidate variables to a regression model, the best subsets approach tests all possible combinations of the predictor variables, and then produces a list of the «best» model(s) according to some statistical criteria (in this case, Akaike information criterion (AIC), which balances the trade-off between the goodness of fit of the model and the simplicity of the model).

Further data analysis included exploration of the relative contribution of various potential drivers (e.g. chemical and climatic variables) in the three time periods (i.e. the full record (1990–2016); the first 15 years (1990–2004; referred to as «early»); and the last 15 years (2002–2012; referred to as «late»)). We used the realimpo package in R [41] to compute the proportion of R² that each predictor variable contributes using the lmg metric [42]. The technique uses variance decomposition to average the improvement in R² gained by adding each candidate predictor to multiple linear regression models of various sizes and 'quantifies the relative contributions of the regressors to the model's total explanatory value by averaging sequential sums of squares over orderings of regressors' [43]. The models were run with six variables (the four primary dynamic drivers, plus two interactions between the chemical variables and median (Ca + Mg)), for a total of 720 model iterations (for each time period). The relative importance of each variable is presented as a share of the unexplained variance, and together sum to 100%. For the purpose of our analysis, we included all possible forms of each dynamic driver, and summed their lmg values to represent their overall effect (e.g. for SO₄, we summed the lmg values for ΔSO₄ and the interaction of ΔSO₄ and median (Ca + Mg)).

We used publicly available estimates of global runoff data (Global Runoff Data Centre; www.bafg.de/GRDC/EN/Home/homepage_node.html; https://geoportal.bafg.de/mapapps/resources/apps/GRDC_FWF/index.html?lang=en) for three European and three North American drainage basins in combination with median DOC concentration values to calculate changes in riverine DOC export to coastal regions. Drainage basins with substantial regional coverage (Baltic Sea, Norwegian Sea, North Sea; Gulf of St Lawrence, North East Shelf, Scotian Shelf) of sites in our database were included, while drainage basins with limited data (Barents Sea and Newfoundland Shelf) were not. Averaged annual runoff data for 1961–2016 were compiled, and runoff was converted to mm yr⁻¹ using the reported area for the drainage basins. None of these regional annual runoff time series showed a long-term trend (figure S4 (available online at stacks.iop.org/ERL/16/104009/mmedia)). Sites were attributed to drainage basins based on location. DOC flux at the start of the period of interest (1990–2006) was estimated as the median DOC concentration (across all sites in the basin) for the initial five years (1990–1994), multiplied by the mean annual runoff; these fluxes are assumed to represent flux values in 1992. Estimated changes in DOC export were calculated by applying the median rate of DOC change for all sites in each basin to the estimates of flux in 1992, and then multiplied with time. Thus, the changes in DOC flux we present are driven entirely by the measured changes in DOC concentrations, with export of water assumed to be unchanged during the period of interest.

Over the period 1990–2016, ΔDOC was overwhelmingly positive (figure 1; supplementary information (SI), figure S1), with 383 of 426 sites showing positive slopes (269 significant, p < 0.05). Only 43 sites showed negative slopes, four of which were significant (p < 0.05). As demonstrated elsewhere [1, 4], the magnitude of ΔDOC depended strongly on the median DOC concentration (figure S2). Long-term average concentrations of DOC can be considered as surrogates for catchment carbon pools [44, 45], while also being partly dependent on the extent of dilution by precipitation [4, 35, 45].

Zoom In Zoom Out Reset image size

Overall, median ΔDOC was higher in the 2002–2016 period than in 1990–2004, but the median proportional rate of change, expressed as %ΔDOC (100 × ΔDOC divided by median DOC concentration), was not (table 1). Thus, across the full dataset, rates of browning have continued at a similar pace to earlier [1]. However, browning rates accelerated in North America whereas they slowed in Europe (figure 2). For the entire period, %ΔDOC in Europe was approximately twice as high as in North America (table 1).

Zoom In Zoom Out Reset image size

The contrasts in regional patterns and rates of browning call for a further analysis of potential controls. We used a best subsets approach to multiple regression to identify critical drivers of ΔDOC (1990–2016), starting with a comprehensive list of chemical and climatic independent variables and their interactions (table S1). Using the AIC we observed that model fitness plateaued at six or seven variables. The best models that could be achieved using fewer variables (e.g. <4) consisted exclusively of chemical variables.

The seven-variable model (six predictors and one interaction) with the lowest AIC explained 67% of the variance in ΔDOC (table 2) and emphasized the importance of median DOC and the chemical drivers (ΔSO₄ and ΔCl) for ΔDOC. The model coefficients imply that declines in SO₄ and Cl deposition (figure 3) exert positive effects on DOC concentrations, in agreement with earlier findings [1]. The interaction between median divalent base cation concentrations and deposition implies that reduced deposition has a greater effect in acid-sensitive catchments (i.e. those with lower base cation availability), consistent with the soil mechanisms outlined above [13]. The model captures most of the regional variation in DOC trends as illustrated in boxplots of region-wise residual variation in ΔDOC (figure 4). One exception is Newfoundland, where high DOC values in the early 1990s are not adequately captured in the model (figure S3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two climate variables representing change in summer precipitation and temperature (Δsummer_P and Δsummer_T) exert small but significant, and opposing, effects, once the chemical drivers have been accounted for. Trends in both variables are predominantly upward (table 1 and figure 5), although more consistently significant for temperature than for precipitation. Increasing Δsummer_P is positively associated with ΔDOC, whereas the opposite is true for increasing Δsummer_T.

Zoom In Zoom Out Reset image size

We tested for evidence of an increased influence of climate on DOC trends, using a variance decomposition method to estimate the relative importance of all candidate predictors in multiple regression models for each period (section 2 and table 1). As we were primarily interested in the influence of dynamic predictors, we forced the analysis to include median DOC as a first step, and then calculated each driver's contribution to the remaining variance (figure 6). For the entire data period, the relative importance analysis reinforces the conclusions from best subsets multiple regression (table 2) that the influence of the chemical drivers exceeds that of temperature and precipitation. However, the proportion of the variance explained by climatic versus chemical drivers changes. In the early period, circa 90% of the explained variance in ΔDOC by these two drivers was attributable to ΔSO₄ and to a lesser extent ΔCl, implying that deposition reductions had driven DOC increases overwhelmingly. In the late period, however, the sum effect of trends in climate variables (predominantly precipitation) accounts for 48% of this variance, i.e. very similar to the sum effect of trends in chemical variables (52%). This change in the relative importance of climate variables is striking and suggests that these ecosystems may be moving towards a situation where climate drivers will dominate spatial variation in size and direction of DOC trends.

Zoom In Zoom Out Reset image size

The accelerating and decelerating rates of browning in North America and Europe, respectively, are highly consistent with differences between regions in the temporal behavior of the dominant drivers (table 1 and figure 5). Most importantly, SO₄ declines have slowed in Europe, but have accelerated (i.e. become more negative) in North America. Additionally, the substantial reductions in Cl concentrations in Europe in the early period did not continue in the late period, while changes in Cl concentrations reversed from slightly positive to negative in North America reflecting both the variable nature of seasalt deposition and changes in the importance of HCl deposition. Summers in North America shifted from a drying trend in the early period to a wetting trend in the late period, while in Europe summer rainfall increased in both periods. Summers became warmer on both continents during the early time period, but this trend continued into the late period in North America only. We conclude that the contrast in acceleration of browning in North America and deceleration in Europe is driven by differential trends in cleaner air, seasalt deposition and summer precipitation.

Using long-term averages of riverine water transport from publicly available data (figure S4) and median concentrations of DOC from our dataset, we estimated how riverine DOC transport to coastal regions from large drainage basins encompassing our study sites might have changed over the course of our study (table 3). This analysis suggests that the total export of DOC from these areas has increased by 27%, i.e. from 8.0 to 10.7 Tg C yr⁻¹ between 1992 and 2016 (table 3). Grouped by continent, the increase was largest in Europe (Europe: +1.7 Tg C yr⁻¹ or +35%; North America: +0.7 Tg C yr⁻¹ or +16%) which agrees with the highest absolute DOC trends and concentrations found here.

Our analysis demonstrates that recent browning rates in headwaters of many European and North American landscapes are comparable to those reported earlier [1], i.e. at rates between 0% and 2% annually. We found no differences in relative DOC trends in the early and late period, and thus no change in browning rates. Most sites present in dataset from Monteith et al [1] were also included here, with the exception of Finland (number of sites declined from 157 to 26) and Ontario (number of sites declined from 32 to 15). Despite considerable differences in vegetation cover, sampling frequency, exposure to air pollution and dominating climate, simple linear models with a limited number of variables could explain two thirds of the observed variation in DOC trends. The relationships we find in this study are consistent with those reported in 2007 [1] with regard to the importance of chemical controls on DOC. However, we detected a stronger effect of climate for the 1990–2004 period, possibly because of the greater spatial resolution of available modeled climate data that was not available in the mid-2000s (i.e. 0.5 × 0.5 vs 5 × 5 degree grids). The larger impact of climate in the 2002–2016 period strongly suggests that climatic factors have become more important for browning, an observation supported elsewhere [4].

While climate change effects on biogeochemical processes controlling carbon cycling are expected to intensify, pollutant S and Cl levels in atmospheric deposition are now so low in many areas in Europe and North America that there is limited potential for further decreases [22]. Seasalt deposition and precipitation patterns appear to be related to atmospheric circulation patterns such as the North Atlantic Oscillation [25], which are inherently complex to forecast [46], rendering predictions of storms that generate seasalt aerosols very uncertain. By contrast, scientific consensus on global warming is beyond dispute [47], with both flooding and droughts expected to become more frequent. Land-use change such as forest expansion [48] could influence DOC export in some areas, but was not found to influence DOC trends in the UK sites in our study [35]. Land use evidently cannot explain either the consistency or the spatial distribution of increases that are observed across an extensive and diverse range of sites included in our dataset, many of which are located in alpine or subarctic areas with little or no forest or have grasslands as the dominant land cover.

The positive correlations we found between precipitation and DOC may be related to the increased proportion of superficial runoff through organic-rich topsoils during periods of high precipitation [49]. Contrary to expectations, we found no evidence that warming has contributed to browning. Instead, we found a small, but significant, damping of DOC increases by higher summer temperatures. Studies of spatial variation have found positive relations between temperature and lake DOC [26], while century-long DOC-reconstructions from lake sediments suggest that higher precipitation and temperature both promote DOC, in addition to reduced acid deposition and land use [50]. The latter study [50] highlights that impacts of such controls are overlapping and their relative importance is difficult to ascertain. Studies of seasonal variation in DOC often show increased DOC during summer and fall, interpreted as a consequence of the temperature dependency of DOC production [26]. However, warming also promotes the mineralization of organic matter to CO₂, and the effect of summer warming could be to switch the balance from DOC production (i.e. partial decomposition of organic matter) towards CO₂ production [3] (i.e. complete decomposition) [51]. In an observational study covering spatial and temporal variation like the one we report here, we assume that drivers exert the same influence in space and in time, i.e. that spatial differences in acidic deposition produce similar differences in DOC behavior that temporal trends in acidic deposition would. With substantial internal correlation between the climatic and chemical drivers, we cannot know for certain whether the observed negative temperature effect is due to a spatial pattern (for example, if higher summer warming is associated with areas that also show larger changes in SO₄) or whether it presents a shift from soil production of DOC to CO₂.

Drought and excess precipitation have a strong influence on net primary production [3, 52], vertical [53] and lateral [54] carbon transport and aquatic carbon cycling, the latter through the positive relationship between photochemical and microbial processing of DOM and water residence time [55]. The influence of climate on vertical and lateral transport of DOC has implications for the overall carbon balance of terrestrial ecosystems. Depending on climatic factors, different proportions of the DOC mobilized from organic soils are likely to be retained in mineral soils, exported to surface waters and re-sequestered in lake and marine sediments, or mineralized in aquatic systems and either incorporated in the ocean inorganic carbon pool or returned directly to the atmosphere. In the context of drought, lower DOC transport in a warming world could be associated with higher soil CO₂ emissions.

Land to sea transport of carbon by rivers is not included in the most recent global carbon budget [56] although its magnitude (i.e. 1.0 Pg C yr⁻¹) is significant compared with the land carbon sink and is likely to have increased significantly since the pre-industrial era [57, 58]. Our estimate suggests an increase of 27% (2.3 Tg C yr⁻¹) from large watersheds in Europe and North America in the course of less than three decades. The average areal export of DOC in our study was between 1.5 and 3.4 g C m⁻² yr⁻¹ which is in the range reported for tundra, taiga and boreal ecosystems [58–60]. Inferring lateral DOC export from annual values of DOC concentrations and freshwater runoff results ignores the often disproportional contribution of snowmelt and precipitation events to DOC export [61], suggesting that our estimate of lateral C export could be too low. However, DOC retention and loss during transport from headwaters to coastal waters by aquatic processing can also be significant [55], suggesting a possible overestimation of our coastal DOC input value. Our assessment also omits any additional DOC sources that may be present in the lower reaches of river systems, for example those associated with agricultural, urban and industrial inputs and does not take interannual variations in DOC and discharge into account. Improved estimations of headwater export and subsequent aquatic processing would be valuable for constraining the role of aquatic ecosystems in the global carbon cycle, in the context where climate is becoming a more important control of DOC trends.