From greening to browning: Catchment vegetation development and reduced S-deposition promote organic carbon load on decadal time scales in Nordic lakes

Increased concentrations of dissolved organic carbon (DOC), often labelled “browning”, is a current trend in northern, particularly boreal, freshwaters. The browning has been attributed to the recent reduction in sulphate (S) deposition during the last 2 to 3 decades. Over the last century, climate and land use change have also caused an increasing trend in vegetation cover (“greening”), and this terrestrially fixed carbon represents another potential source for export of organic carbon to lakes and rivers. The impact of this greening on the observed browning of lakes and rivers on decadal time scales remains poorly investigated, however. Here, we explore time-series both on water chemistry and catchment vegetation cover (using NDVI as proxy) from 70 Norwegian lakes and catchments over a 30-year period. We show that the increase in terrestrial vegetation as well as temperature and runoff significantly adds to the reduced SO4-deposition as a driver of freshwater DOC concentration. Over extended periods (centuries), climate mediated changes in vegetation cover may cause major browning of northern surface waters, with severe impact on ecosystem productivity and functioning.


B) Additional sampling site and source data description Lake chemistry data
The data in this analysis are obtained from a Norwegian lake monitoring program from 70 lakes covering the Norway mainland and with annual samples from 1986 to 2013. These lakes constitutes a subset of the 1000-lake acidification survey of 1995 (Henriksen et al. 1998, Ambio 27:80-91), and represents acid sensitive, headwater lakes on granitic or gneissic bedrock with negligible local pollution sources. Water samples are collected annually at the outlet after the autumn circulation period and where analysed at the Norwegian Institute of Water Research (Skjelkvåle et al. 2001,Hydrol. Earth Syst. Sci. 5:327-338). The original dataset (available from URL http://vannmiljo.miljodirektoratet.no/)contains information from 78 lakes, but eight lakes where excluded from the present study due to incomplete catchment data. Catchments crossing national borders did not have complete data on runoff (provided by The Norwegian Water Resources and Energy Directorate). Table S1: List of study sites (lakes) included in the present study. The table is given with lakeID corresponding to the official lake numbering given by the Norwegian Water Resources and Energy Directorate, http://www.nve.no/en/, latitude and longitude (EPSG:4326) of lake centroid, lake area (square kilometers), catchment area (square kilometres) and mean lake TOC concentration.  . The NDVI is a radiometric measurement of the fraction of photosynthetically active radiation (~ 400 to 700 nm) absorbed by chlorophyll in the leaves of a vegetation canopy (Myneni et al. 1995, Trans. Geosci. Rem. Sens. 33:481-486). Although care should be taken before interpreting NDVI changes directly as net changes in primary production, the index has previously proven a good surrogate of photosynthetic activity and for capturing long-term changes in vegetation (Myneni et al. 1997, Nature 386:698-702;Forkel et al. 2013, Remote Sensing 5:2113-2144. The NDVI3g data set has a pixel resolution of 8 x 8 km. The maximum NDVI value over a 15-day period is used to represent each 15-day interval to minimize bias due to cloud contamination and variations in atmospheric turbidity, scan angle, and solar zenith angle (Holben 1986, Int. J. Rem. Sens. 7:1417-1434Pinzon and Tucker 2014, Remote Sensing 6:6929-6960). This scheme results in two maximum-value NDVI composites per month. Here, we used the mean NDVI3g record for the main snow free period (June -August) averaged annually from 1982 to 2011 and over each catchment area. Other measurements, such as maximum yearly NDVI3g and June -August maxima gave similar results in terms of temporal dynamics, and were therefore not included in further analyses. Time series of S deposition at the catchment level was interpolated from EMEP deposition raster temporal composite (1980,1985,1990, and annually from 1995 to 2011) (Schulz et al. 2013, Transboundary acidification, eutrophication and ground level ozone in Europe in 2011. Norwegian Meteorological Institute) using recalculations if available. Total S deposition was used in the present study, comprising both dry and wet deposition. Each type of deposition was modelled using regression kriging using EMEP precipitation as a linear regression model predictor, with residuals fitted using a linear variogram model with no nugget. Using the variogram model and linear model, deposition flux over the catchment area was predicted using block kriging with the target geometry being the catchment polygon. Missing values from 1985-1990 and 1990-1995 where estimated within each catchment using loess regression.

C) Extended output and additional analyses Data standardization and software
All pre-processing and analyses was done in R v. 3.2.1 (R Core Team 2015). Prior to the analyses, TOC where log transformed, and the data where standardized by detracting mean and dividing on standard deviation using the scale function of the base library.
We first tested for overall temporal trends in the time series of lake TOC and each of the catchment drivers using a regional Kendall trend test using the rkt library (Marchetto 2015, rkt: Mann-Kendall Test, Seasonal and Regional Kendall Tests). Each individual lake was used as block (Helsel and Frans 2006, Env. Sci. Tech. 40:4066-4073).
To test for effects of vegetation development ( New York). This resulted in a full model including lake-ID as random intercept and year nested within lake-ID as random slope (???AIC > -76.65 in favor of the selected random structure). By introducing lake ID as random factor, we model between-lake variation in TOC caused by catchment variables not considered in the current study. For example static land cover properties that not are expected to change on the investigated time frame such as bogs and wetlands, catchment size, lake size. Temporal autocorrelation (AR1) was initially also tried, by entering year nested within lake-ID, but proved redundant in the initial selection of the random structure. There was no signs of multicollinearity (variance inflation factors, all VIF < 1.96, (Zuur et al. 2009) and maximum correlation between predictors where -0.61.
Catchment vegetation is expected to influence lake TOC with time lags (see main text). Due to catchment-specific properties such as fractions of alpine areas, forests and bogs, as well as inter-annual variability in climate and runoff, the time-lag between terrestrially fixed C (inferred as NDVI) and export of TOC may vary substantially. Based on visual inspection of cross-correlation plots between the dependent variable (TOC) and time lags of the explanatory variables, it was apparent that NDVI and runoff affected lake TOC with a considerable time lag. In comparison, there were no clear time-lags in the effect of S deposition or temperature.
Figure S1: Inspecting for delayed response (time-lags) between drivers and response Here, we do a visual inspection of cross-correlation plots between the dependent variable (TOC) and time lags of the explanatory variables. The correlations are given as means across all lakes, and solid line represent a loess fitted curve for illustrational purposes.
NDVI and runoff affected lake TOC a considerable time lag. In comparison, there were no clear time-lags in the effect of Sulphur deposition or temperature. In order to avoid data-dredging issues, we selected a priori a one year time lag for NDVI and runoff as the input for our analyses. We subsequently performed sensitivity analyses on the effect of time-lag chosen by re-running the final model with different time-lags (see below). The main predictions were mainly insensitive to choice of time-lag, and hence only 1 year time lags are presented in the results of the main text. S deposition and temperature where entered without time lags.

Model selection and multi-model inference
Model selection and model averaging on fixed effect structure where done by model comparison using the MuMIn library (Grueber et al. 2011, J. Evol. Biol. 24:699-711).
We compared all possible combinations of fixed factors and ranked candidate models according to the Akaike information criterion (AIC), and also calculated delta AIC ( i) in relation to the highest ranked model, as well the Akaike Weigths (wi). There was no clear top-ranked candidate model, and we applied Akaike weigth-based averaging over the 95% confidence model set (cumulative AIC weights of models 0.95) for estimating coefficients for the candidate models as well as their 95% confidence intervals. The relative influence (RI) of each variable was given as the summation of wi across all models including that variable in the 95% confidence model set (Johnson and Omland 2004, Trend. Ecol. Evol. 19:101-108).

Rerunning model exluding outliers
Outliers identified in residual plots (see above, standarized lnTOC < -2.5). The outliers are clusters of low TOC values. We have no knowlegde about potential causes for the outliers, and hence no a priori reason for exluding them from the analyses. One potential hypoteshis is, however, that this results from TOC beeing close to the detection limits of the instruments. We consecvently rerun the model selection process without these datapoints to check for influence. Although the relative importance of NDVI became sligtly lower, the main results of the model selection is similar, albeith a sligthly lower relative importance of NDVI. Hence, our conclution is that these outliers don't have any significant effect of our overall conclutions.

Re-running model selection with different time lags for NDVI and runoff(Q)
We performed sensitivity analyses on the effect of chosing of time-lag by re-running the final model with different time-lags. Model selection was run through all combinations of NDVI and runoff(Q) lags from 0 to 5, and the relative importance of NDVI in the model recorded. The main predictions were mainly insensitive to choice of time-lag, and hence only 1 year time lags are presented in the results of the main text. S deposition and temperature where entered without time lags. Figure S4: Visual presentation of support for NDVI effect using different time lags for NDVI and runoff. The effect of different time-lags on the relative importance of NDVI visualized by plotting the relative importance of NDVI along different time lags of runoff (Q) and NDVI. Range in relative importance is from 0.22 (dark green) to 0.99 (red). All time-lag combinations investigated included NDVI in the top ranked models (delta AIC < 2).

Substituting S deposition with SO4
We then illustrate the effect of substituting SO4 with NDVI in the original model. No major differences on model selection results by subtituting S deposition with SO4 measured in water. NDVI had for all practical purposes the same importance as in the orignal model using S deposition model. Hence, the two explanatory variables (S deposition and SO4) yielded comparable results with respect to model output.

D) Trend analyses
In this part of the supplementary, we investigated for the relative correspondence between the temporal trends in lake TOC and catchment drivers by using the Theil-Sen's slopes for lake TOC as dependent variable and slopes for mean catchment specific NDVI, runoff, atmospheric S deposition as explanatory variables using generalized least-squares (gls). We also tested for the inclusion of spatial autocorrelation structure comparing models with and without (fitted by REML) using Akaike's information criterion (AIC) (Zuur et al. 2009). However, no support for inclusion of a spatial autocorrelation was apparent (∆AIC > 4.00 in support of model without spatial autocorrelation). The the tests is also repeated using relative trends (% change), as well as absolute change for sensitivity puroposes. We additionally present analyses where S deposition where substituted with SO4 as dependent variable, and finally re-run with trends expressed as percent change instead of slopes.  Using trends (Theil-Sen's slopes) but with measured SO4 in lake water samples as predictor   Using relative slopes (% change) with S deposition as driver Figure S6: Site-specific percent change of S deposition, runoff, summer NDVI, and lake TOC concentration, relative to the mean of all years available at each lake and catchment. Figure