Can the heterogeneity in stream dissolved organic carbon be explained by contributing landscape elements?

. The controls on stream dissolved organic carbon (DOC) concentrations were investigated in a 68 km 2 catchment by applying a landscape-mixing model to test if downstream concentrations could be predicted from contributing landscape elements. The landscape-mixing model repro-duced the DOC concentration well throughout the stream network during times of high and intermediate discharge. The landscape-mixing model approach is conceptually simple and easy to apply, requiring relatively few ﬁeld measurements and minimal parameterisation. Our interpretation is that the higher degree of hydrological connectivity during high ﬂows, combined with shorter stream residence times, increased the predictive power of this whole watershed-based mixing model. The model was also useful for providing a baseline for residual analysis, which highlighted areas for further conceptual model development. The residual analysis indicated areas of the stream network that were not well represented by simple mixing of headwaters, as well as ﬂow conditions during which simple mixing based on headwater watershed characteristics did not apply. Speciﬁcally, we found that during periods of baseﬂow the larger valley streams had much lower DOC concentrations than would be predicted by simple mixing. Longer stream residence times during baseﬂow and changing hydrological ﬂow paths were suggested as potential reasons for this pattern. This study highlights how a simple landscape-mixing model can be used for predictions as well as providing a baseline for residual analysis, which suggest potential mechanisms to be further explored using more focused ﬁeld and process-based modelling studies.


Introduction
Dissolved organic carbon (DOC) is a key constituent in surface waters as it has fundamental implications for the ecology and biogeochemistry of aquatic ecosystems. The important role of stream DOC has resulted in several recent investigations to better understand the mechanisms of DOC regulation across temporal and spatial scales (Tank et al., 2012;Temnerud and Bishop, 2005). A general finding has been that the variability of stream DOC concentrations within and between adjacent streams can be as large as the variability found on a regional or even global scale . Although much of this variability can be explained by the occurrence of organic soils in the catchments (Creed et al., 2003;Walker et al., 2012), peatlands alone do not explain the large spatial heterogeneity of DOC in the landscape (Mattsson et al., 2009;Ågren et al., 2007).
The highest surface water DOC concentrations in midto high latitudes can be found in the boreal biome (Mulholland, 2002). Simplified, this regional pattern occurs because of climatic and hydrological conditions including relatively low temperature and high soil moisture, limiting the mineralisation of litter/soil organic carbon relative to plant growth and litter production. This leads to a buildup of organic matter which is subsequently available for export to streams as DOC . Despite these generally high DOC concentrations in the boreal landscape, there is a large spatial heterogeneity in stream water DOC derived from different landscapes that varies because of catchment characteristics, topography and size (Temnerud et al., 2007; Published by Copernicus Publications on behalf of the European Geosciences Union.  Tank et al., 2012). While the accumulation of organic matter in unforested mires makes them the major source of DOC in the boreal landscape (Rantakari et al., 2010;Ågren et al., 2007), forested areas, which generally have the greatest areal extent in the boreal biome, also contribute large DOC concentrations because of the presence of organic-rich riparian soils Knorr, 2013).
The relative proportion of mires and forest in the landscape can be used as a first-order approximation to predict the stream DOC concentration in small streams (Aitkenhead et al., 1999;Laudon et al., 2012). However, as the catchment size increases from headwaters to meso-scale catchments, so does the complexity of the contributing factors controlling stream water chemistry (Bloschl and Sivapalan, 1995). This increased complexity can be related to new/different contributing landscape features becoming increasingly common downstream at lower elevations, but also because there may be scale-dependent processes that can have considerable effects on the stream DOC concentrations as the rivers grow, for example changing flow paths (Cey et al., 1998) or effects of landscape structure (Pacific et al., 2010).
Another characteristic feature of DOC is the large temporal variability related to hydrological events, seasonal differences and inter-annual conditions (Dawson et al., 2011). Hydrology has a first-order control on DOC concentrations in individual catchments (Hinton et al., 1997;Laudon et al., 2011;Raymond and Saiers, 2010). Drying and re-wetting of catchment soils (Köhler et al., 2008), soil temperature (D'Amore et al., 2010), winter climatic conditions  and antecedent conditions controlling the pool of sorbed, potentially soluble organic carbon Yurova et al., 2008) can affect DOC concentrations on an event, seasonal and annual timescale. This temporal variability adds to the spatial complexity of DOC concentrations, as different catchment characteristics can differ in response to hydrological and climatic forcing depending on catchment soils, vegetation and topography. Furthermore, depending on the spatial configuration of the landscape, the residence time of water in the surface water network can moderate or exaggerate the response in downstream locations in ways that are not easily predictable.
Because of the large complexity of factors controlling stream DOC concentrations we tested a simple conceptual landscape-mixing model as a predictive and diagnostic tool on a large nested boreal stream data set, to better understand how DOC is regulated during different seasons and across scales. The main objectives of this study were (1) to test if the spatial heterogeneity of stream DOC concentrations can be explained by the major contributing landscape elements; and (2) to use a residual analysis as a diagnostic tool and a learning framework for further development of our conceptual understanding. To answer these questions we used a landscape-mixing model on the 68 km 2 boreal Krycklan catchment. A landscape-mixing model (Cooper et al., 2004;Cooper et al., 2000;Evans et al., 2001) offers a sim-ple approach to modelling stream biogeochemistry by lumping processes into dominating landscape elements that can be used to examine if the DOC concentration is simply due to the conservative mixing of contributing sources. Firstly, we investigate how the landscape-mixing could be used to predict DOC, using several model assessment criteria. Secondly, we analyse the model residuals to investigate model performance and thereby answer the question of where in the landscape simple mixing of stream water does not adequately characterise stream DOC behaviour. By running and validating the model on data sampled on seven different occasions, we also addressed the question of whether the landscapemixing model performed better at certain times of the year, or under certain flow conditions. Using this simple conceptual approach the ultimate goal was to provide new insights into the mechanisms regulating stream DOC and how these may vary across a landscape and during different times of the year.

Study catchment
The 68 km 2 Krycklan catchment (64 • 16 N, 19 • 46 E) was used as a study catchment for modelling the spatial variability of DOC in the stream network . The Kryklan catchment is a glaciated forested catchment (forest cover 87 % and peatland cover 9 %). The forest is dominated by Scots pine (Pinus sylvestris) and Norway spruce (Picea abies) with an understorey dominated by ericaceous shrubs, mostly bilberry (Vaccinium myrtillus) and lingonberry (Vaccinium vitis-idaea) on moss mats of Hylocomium splendens and Pleurozium schreberi (Forsum et al., 2008). Quaternary deposits of till, peat and fine sorted sediments are the dominant overburden (Fig. 1). The peatlands are classified as forested (2/3) or open (1/3). The peat is dominated by Sphagnum species and consists of mostly minerogenic, acid and oligotrophic mires with varying proportions of microtopographic units (e.g. strings/lawns). Rock outcrops and thin soils are common on hilltops. The region is affected by isostatic postglacial rebound. The highest postglacial coastline crosses the catchment at approximately 256 m a.s.l., and 45 % of the catchment is situated above the former coastline. In the lower lying parts of the catchment, there is an old postglacial delta with deposits of mostly silty sediment. While there are a number of small lakes in the catchment, the overall lake coverage is small at 0.6 %. Human impact is low and agricultural land covers only 2 % of the catchment. The bedrock consists of 94 % sedimentary rocks (Precambrian metagreywacke) with smaller patches of basic volcanic rocks and acid volcanic rocks, covering 3 and 4 % respectively. The site conditions are characterised by long winters, with snow cover typically from November to the beginning of May. The 30-year (1981-2010)  1.8 • C and the annual precipitation is 640 mm, of which approximately half enters streams as runoff (Oni et al., 2013).

Stream water sampling
Stream water was sampled at 115 sites throughout the catchment on seven occasions from May 2003 to Sept 2008 during different seasons and hydrological conditions (Table 1 and Fig. 2). The sampling campaigns were designed to take a "snapshot" of the spatial variability of the stream network, and on each occasion, sites were sampled during a single day, except during winter baseflow, where sampling extended over a week. While 115 different site locations were sampled in total, the number of sites sampled in any particular survey varied between 73 to 89. A subset of 42 sites was sampled on all seven occasions.
Discharge was measured in a second-order stream in the central area of the catchment (called Svartberget or C7 in previous studies; Laudon et al., 2007) at a 90 • V-notch weir located inside heated housing. Pressure transducers connected to Campbell scientific data loggers, USA or duplicate WT-HR capacitive water stage loggers, Trutrack Inc., New Zealand were used to record the water level. Using established rating curves the water stage was used to calculate discharge. We make the simplifying assumption that the specific discharge is the same throughout the catchment. The uncertainty this assumption introduces has been calculated to be on average at most 12 % (Ågren et al., 2007), but can be higher under particularly low flow conditions . The water samples were collected in acid-washed and sample-rinsed high density polyethylene (HDPE) bottles (Embalator Mellerudplast, Mellerud, Sweden) and were stored frozen until they were analysed for DOC using a Shimadzu TOC-VCPH/CPN analyser (Shimadzu, Kyoto, Japan).

Watershed characteristics
Lidar (Light detection and Ranging) measurements of the catchment have been made at a point density of 3.3-10.2 measurements per m 2 . These data were used to generate a 0.5 m high-resolution DEM. For hydrological modelling the DEM was aggregated to a 5 m resolution. In order to make the DEM flow compatible it was manually corrected where bridges and road culverts obstructed the flow algorithm, and all sinks were filled. The catchment delineation was then derived automatically from the DEM using ArcGIS 10.0. Care was taken to ensure that the catchment delineation was correct for all 115 catchments, and manual adjustments were made to the DEM in questionable sections based on a 3-D version of the 0.5 m DEM combined with field observations. For each sub-catchment the catchment characteristics were derived using map data. DOC was modelled from the surficial geology cover based on the Quaternary deposits map (1 : 100 000) (Geological Survey of Sweden, Uppsala, Sweden). Additional catchment characteristics were derived for all sub-catchments for potential use as covariates in the residual analysis. These characteristics included stream order, catchment area, slope, topographic wetness index (TWI MD8 ) , proportion above the highest coastline, as well as the land cover from the "road map" (1 : 100 000) and the "property map" (1 : 12 500) (Lantmäteriet, Gävle, Sweden). With the aid of IR orthophotos, combined with a detailed forest inventory, the whole catchment has been divided into 1751 forest stands (areas with a similar mix of known tree species and age). Based on the Table 1. Discharge (mm day −1 ) for each sampling occasion and the estimated end-member concentrations (mg L −1 ) from bootstrapping (n = 15). R 2 is the R 2 from the bootstrapping procedure (Eq. 1). To the right is the uncertainty in modelled concentrations expressed as the coefficient of variation (CV) from the Monte Carlo analysis. Lidar measurements and regression models with field observations, detailed maps were constructed providing, for each 10*10 m pixel, forest stand height, birch (Betula spp) biomass, lodgepole pine (Pinus contorta) biomass, Norway spruce biomass, Scots pine biomass, total biomass and mean forest stand age. Averages of all the forest variables were calculated for each sub-catchment.

End members and landscape-mixing modelling
The landscape-mixing model, which was based on Cooper et al. (2004) and Cooper et al. (2000), predicts water chemistry throughout a stream network from landscape properties. The model is based on the assumption that the variability within a landscape type is smaller than between landscape types, and that different landscape elements generate different solute concentrations. These landscape concentrations are estimated from sampling data at stream locations draining subcatchments with known upstream proportions of each landscape type. A detailed DEM (digital elevation model) with 5 m resolution and the presence of many sampling sites in our study allows us to work with the actual sub-catchments and at a high resolution. We used a statistical approach to calculate the end-member concentrations from the different landscapes. In order to more easily compare model performance between all seven sampling occasions we selected headwater catchments that were sampled on all occasions as the data set for model parameterisation. Small catchments with more uniform landscape characteristics will tend to have concentrations which are closer to the different end members and more representative of sources, while larger catchments, through mixing of upstream sources, show a reduced variability (Temnerud and Bishop, 2005). We therefore selected only catchments with area less than or equal to 3 km 2 for model parameterisation. Fifteen catchments fulfilled both criteria (sampled on all occasions, size ≤ 3 km 2 ); the remaining sampling sites were used to assess model performance, particularly to test the simple mixing hypothesis.
Previous research in the catchments has identified three landscape types which are expected to give rise to contrasting stream water chemistry, including DOC concentrations Buffam et al., 2008). We have termed these landscape types "peat", "till", and "sorted sediments", based on the corresponding surficial geology deposits underlying each landscape. The variation in surficial geology influences other landscape characteristics including weathering rates and drainage, which in turn influence soil formation, vegetation, subsurface hydrologic flow paths and rates, and riparian zone formation. All of these are expected to influence DOC, thus the surficial geology categories serve as a useful tool for categorisation. From the surficial (Quaternary) geology map each of the 115 catchments was classified by relative proportion of (A) peat, (B) till (this also includes the "thin soil" class which in essence is a shallow layer of till on bedrock), (C) sorted sediments (silt, sand and glaciofluvial alluvium) or (D) "other" (lakes and rock outcrops). Based on the 15 selected headwater sites in the construction data set, a regression model (Eq. 1) was constructed to calculate the end-member concentrations for each landscape type and on each sample date by multiplying the concentration with the areal coverage for each landscape type (A-D). By setting the intercept to 0 in the model and using the areal coverage of the landscape types in proportions (0-1) instead of percentages, the estimates (A-D) were expressed directly as the end-member concentration for DOC in mg L −1 for each landscape type.
Because of non-normal distribution of data, a bootstrapping approach was used to solve the "landscape concentrations" iteratively. The bootstrapping procedure was done by sampling with replacement to generate samples of the same size as the original data set. Under this procedure, a random number of streams were deleted from the data set, from the remaining data set some streams were then included twice, or more, until the data set again comprised 15 streams. Slopes and constants were calculated for every new data set, then the randomisation process was repeated 1000 times. Finally, mean concentrations were computed for each landscape component and used in the model (Table 1). This method has the additional benefit that it provides an estimate of the uncertainty in the end-member concentrations. Based on the repeated runs, the standard errors, confidence intervals, and correlations were calculated for each end-member concentration. The uncertainty in the calculated end-member concentrations was later used to analyse the total uncertainty of the models. All bootstrapping calculations were done in PASW Statistics 18 (SPSS Inc.). Initially, the bootstrapping procedure sometimes generated unrealistic estimates. To overcome this, constraints were set on the endmember concentrations. Soil water data from the catchment were used as constraints for concentrations of each landscape type. For peat, lower and upper limits of 4 and 84 mg L −1 were set, based on measurements from groundwater wells in a wetland in the catchment (Yurova et al., 2008). For till the acceptable range was set to 1-97 mg L −1 given the variability in lysimeter measurements from 10 soil profiles in till-derived soils in the catchment . In fine sorted sediments the constraint was set to 1-46 mg L −1 given the variability in lysimeter measurements from three soil profiles in the fine sorted sediment-derived soils in the catchment . In the first attempt, the endmember concentration for landscape type "other", consisting of lakes and bare rock (D in Eq. 1), was calculated. The evaluation showed that the end-member concentration for D was extremely variable and uncertain and including these values did not improve the fit for the overall model. Because of this uncertainty and since class "other" had such a minor areal coverage (on average about 2 % and at maximum below 10 % coverage; Fig. 3), the parameter D in Eq. (1) was set to 0.

Landscape-mixing modelling in GIS
The high-resolution DEM facilitated modelling of DOC concentrations every 5 m throughout the entire stream network using the landscape-mixing model and ArcMap 10 hydrological modelling tools. Using a weighting raster containing the end-member concentrations for the aggregated surficial geology map (aggregated into the four classes) when performing the flow accumulation calculation, the DOC export from each cell was calculated. The DOC export was then divided by estimated discharge to calculate the DOC concentrations for all 5 × 5 m cells in the landscape. The modelled DOC concentration for the sampling sites could then be extracted. The modelled DOC values were compared to the measured values for the respective site on each sampling occasion. A layer showing the modelled DOC concentrations for every 5 m section of the stream network could also be displayed (Fig. 1).

Model validation
Model performance was assessed using data from the sites that were not used for model construction. We calculated several measures (Table 2 and Fig. 4). Root mean square error (RMSE) has the benefit that it gives the error in units of mg L −1 . To standardise RMSE we calculated the RMSEobservation standard variation ratio (RSR). A low RSR indicates a better model and values below 0.7 are considered a satisfactory model (Moriasi et al., 2007). As a measure of the average tendency of the modelled values to be larger or smaller than observed values, the percent bias (PBIAS) was calculated. For PBIAS the optimum value is 0, negative values indicate a model overestimation bias and positive values an underestimation of modelled values. We also plotted the measured and modelled values (Fig. 4) and used standard regression measures of R 2 and slope. A slope near 1 indicates that the model is close to the 1 : 1 line, a large diversion from 1 indicates a systematic error in the model. As an example, a slope below 1 means that high DOC concentrations are underestimated and low values are overestimated. The R 2 value indicates the strength of the relationship is between the measured and predicted values, but does not take into account any systematic errors in the slope of the relationship. Nash-Sutcliffe efficiency (NSE) indicates how well the scatter fits the 1 : 1 line; the value of NSE is similar to R 2 in that a value close to 1 indicates a good fit and a value close to 0 indicates a poor fit.

Uncertainty
One source of uncertainty in the model is the representativeness of the 15 selected catchments. To test how this affected the modelled DOC, the bootstrapping routine was rerun using all available sites to calculate the end-member concentrations for the entire data set. The landscape-mixing model was then rerun on new estimates and an evaluation on how that affected RMSE, RSR, PBIAS and NSE was calculated (Table 3). A second source of uncertainty was related to the endmember concentrations. However, by using the bootstrapping method this uncertainty was calculated (Fig. 5). A Monte Carlo analysis was performed to propagate the uncertainties of the end-member concentrations to calculate an overall uncertainty for the modelled values. The overall uncertainty was calculated using 10 000 realisations with random parameters assuming that the uncertainty in A, B and C was normally distributed. The uncertainty was expressed as a coefficient of variation (standard deviation of the modelled values/average of the modelled values) (Table 1).

Residual analysis
The residuals (modelled -measured) were analysed using a multivariate statistical approach, partial least squares   projections to latent structures (PLS), to identify where and when the model failed to reproduce the measured data well. PLS is a method for relating two data matrices, X and Y, to each other by a linear multivariate model (Eriksson et al., 2006b). PLS is similar to principal component analysis (PCA), but instead of extracting the principal components so that they maximise the variance in the X matrix (as in PCA) the PLS method extracts the principal components so that they maximise the correlation between the X matrix and the Y matrix. The strength of the PLS method is the ability to analyse data with "many, noisy, collinear, and even incomplete variables in both X and Y" (Eriksson et al., 2006b, a). The PLS analysis was conducted using the multivariate statistical program SIMCA-P + 12.0.1, Umetrics, Umeå. The first step was to get an overview of the relationship between the response and explanatory variables and the residuals. To achieve this, the residuals (modelled DOC -measured DOC) for all occasions (Y matrix) were related to the catchment characteristics (X matrix) using a single PLS model (Fig. 7). When it was found that the behaviour of the residuals varied according to hydrological conditions, two new refined models were constructed, one for high/intermediate flow (Fig. 8a) and one for baseflow (Fig. 8b). In order to facilitate the interpretation of the graphs in Fig. 8, the PLS models were refined to find the best predictor variables, based on the conditions that the variable coefficient should be significant (95 % confidence interval) and the variable importance of the projection (VIP) should be high (> 1). In SIMCA-P+, for every model, the program also calculates the variable influence for each variable, called variable importance in the projection (VIP). VIP is the sum over all model dimensions. Variables with large VIP, larger than 1, are the most relevant for explaining Y.

Results
The bootstrapping estimates of the end members show that peat has the highest DOC concentrations followed by till and lastly, fine sorted sediments (Table 1 and Fig. 5). Plotting the end-member concentration as a function of discharge for the sampling occasions revealed that the DOC concentrations increased with discharge (Fig. 5). For silt and till this increase could be approximated by a linear relationship. For peat, the curve estimation procedure in PASW suggested a sigmoid curve (p < 0.1). The standard error for the endmember concentration was low for till (on average 2 mg L −1 ) but higher for peat and fine sorted sediments (on average 9 mg L −1 ), where sediment has the relatively highest standard error (Fig. 5). Figure 1 shows an example, from Sept 2008, of the modelled DOC concentrations using the landscape-mixing model combined with GIS and a high-resolution DEM. This shows the strength of this approach during a time when the model performed well and could be used for prediction. With this approach, DOC concentrations can be modelled throughout an entire stream network based on a few headwater observations. It is clear that many of the streams originate in peatlands and have high concentrations initially. As the streams run into the area dominated by till and thin soils the concentrations begin to decrease due to the intermediate concentrations from that composite landscape type. The streams draining the sedimentary area in the valley have the lowest DOC concentrations and as these small streams mix into the larger main stream the concentrations continue to decrease towards the outlet.
The many measures used for evaluating model performance showed somewhat different results ( Table 2). The root mean square error (RMSE) ranged from 2 to 5 mg L −1 . Following the guidelines from Moriasi et al. (2007) the RMSEobservation standard variation ratio (RSR) indicated that only two occasions are considered to be modelled satisfactorily (RSR < 0.70). The mostly negative PBIAS values found (Table 2) indicate a general model overestimation bias. However, using this measure all models performed reasonably well, except for the February 2005 data. The plotted measured and modelled values (Fig. 4) and the slope indicate a systematic bias on all occasions, demonstrating that high DOC concentrations were underestimated and low concentrations were overestimated in the model. The severity of this phenomenon varied and on three occasions the slope was judged to be good, while the other four occasions were judged to be unsatisfactory (baseflow and April 2004) (Table 2 and Fig. 4). According to the model evaluation guidelines by Moriasi et al. (2007), based on the Nash-Sutcliffe efficiency (NSE) measure only two models would be classified as satisfactory (NSE > 0.5). To summarise, the many measures of model efficiency gave different results and contained different information. The most suitable model fit measure depends on the question we are trying to answer. We believe that the RSR (standardised RMSE) and NSE (Nash-Sutcliffe efficiency) give the best overall description of the model performance. Taking into account all model performance measures, the interpretation is that two models performed well (May 2003 andSeptember 2008), one model performed unsatisfactorily (February 2005) and the rest performed satisfactorily.

Uncertainty
Overall, the representativeness of the 15 selected sites was good, but a few outliers were found in the validation data set (Figs. 3 and 4). There were for example a few sites with higher peat coverage than any of the sites that were used for constructing the model (these high peat sites are highlighted as unfilled circles in Fig. 4). To evaluate how this affected the end-member concentrations, all sites (n = 73-89) were used to calculate the end members in Eq. (1). The performance of these new models was then compared to the initial models created from n = 15 sites. For most occasions, the model performance did not change substantially; including all sites to calculate the end-member concentrations only improved the model's predictions by between 5 and 10 % (judged from improvement in NSE) (Table 3). However, in the worst case (February 2005) the improvement was 73 %, indicating that the original construction data set sites were not representative for this occasion. This shows that there is room for model improvement by increasing the number of observations used to calibrate the regression model. However, for this study we wanted to leave as many sites as possible for validation and analyses of the residuals.
The second source of uncertainty was related to the uncertainty of the estimates of the end-member concentrations. Using the uncertainty from the bootstrapping estimates, a Monte Carlo analysis was performed to propagate the uncertainties into an overall uncertainty for the modelled values. As expected, when there were difficulties in constructing a good bootstrapping model, indicated by low R 2 for the model (Table 1, Fig. 6), the uncertainties of the modelled values were high. The model uncertainty also affected its performance (tables 1 and 2 and Fig. 6). February 2005 had the highest uncertainty and had the worst model fit, while the models that performed best (May 2003 and September 2008) had a lower model uncertainty.

Residual analysis
The first PLS model that gave an overview of the data had two significant components (R2Y = 0.35, R2X = 0.57, Q2 = 0.21) (Fig. 7). R2Y and R2X are goodness of fit measures. That means that 57% of the variability in X was used to explain 35 % of the variability in Y . Q2 is an estimate of the predictability of the model. It is calculated by crossvalidation and resembles R 2 in regression models where 0 is poor and 1 indicates optimal predictability. In a PLS loading plot, variables that lie close together co-vary, so the PLS analysis of the residuals (Fig. 7) showed that the residuals clustered based on the discharge of the sampling occasion. The residuals from the high/intermediate flow situations clustered along the first component (black squares in Fig. 7), while the residuals from baseflow measurements clustered higher along the second axis (black triangles in Fig. 7). In order to interpret which variables correlate with high residuals, two different models had to be constructed, one for high/intermediate flow, and one for baseflow conditions (Fig. 8).
Both the model for the high/intermediate flow and the one for baseflow gave PLS models with one significant component. The models were refined to find the best predictor variables, based on the conditions that the variable coefficient should be significant (95 % confidence interval, i.e. significance level of 0.05) and the variable importance of the projection (VIP) should be high (> 1). The models were then rerun on the selected variables to create two refined models ( Fig. 8a and b). The PLS analysis during high/intermediate discharge created a model with R2Y = 0.34, R2X = 0.66, and Q2 = 0.29 (Fig. 8a). That means that the variability in the 23 X variables could be reduced to one component, explaining 34 % of the variability in Y . The nine significant X variables with a high weight explained 66 % of the variability in the extracted component. The PLS loading plot shows the Y weights (c) and the X weights (w * ). The PLS easily handles many covariate variables (Fig. 8a); all X weights that correlate positively to Y weights are different measures related to prevalence of peatlands, and all X weights that correlate negatively to Y weights are different measures of the prevalence of forests. The interpretation of Fig. 8a is that during high/intermediate discharge the landscape-mixing model overestimated the DOC concentrations from sub-catchments with a high coverage of peatlands, and underestimated DOC values in sub-catchments with a high proportion of forest. In contrast, the PLS loading plot for the residuals during baseflow (1 significant component, R2Y = 0.33, R2X = 0.50, Q2 = 0.26) (Fig. 8b) showed that the overestimated DOC concentrations were found in large low-elevation sub-catchments while the underestimated values were those in sub-catchments dominated by peatlands.

Selection of end members for the mixing model
We found that peatlands were associated with the highest DOC concentrations, followed by till and sorted sediments   Fig. 5). Dissolved organic carbon in northern temperate and boreal streams is mostly of terrestrial origin, and peat-containing wetlands are often the major source of DOC (Creed et al., 2008;Dillon and Molot, 1997;Evans et al., 2007;Gergel et al., 1999;Walker et al., 2012). Streams draining the silty sediment area had the lowest concentration, which can be explained by a combination of factors. In Krycklan, catchments underlain by silty sediments are located in the valley bottom of the lower-elevation larger catchments (Fig. 1). A combination of longer flow paths and a high subsurface water transit time can increase the decomposition of DOC (Wolock et al., 1997;Laudon et al., 2007). In addition, a high specific surface area of the fine sorted sediments can lead to increased adsorption to mineral surfaces (Kalbitz et al., 2003).
As previously described, the mixing model is based on the assumption that the variability within a landscape type is smaller than between landscape types, and that different landscape elements generate different solute concentrations. This assumption was true during the high-flow situations, indicated by the separation of the error bars in Fig. 5, but during baseflow there was some overlap of the variability in the end members. Previous research has found that hydrology has a first-order control on the temporal variability of the DOC concentrations in streams (Hornberger et al., 1994;Seibert et al., 2009). Plotting the end-member concentrations as a function of discharge gave a slight positive relationship between DOC and discharge for all landscape types (Fig. 5). We expected, but did not find, a negative relationship between DOC concentration and discharge in mire-dominated catch-ments as suggested by other studies in the study area (Ågren et al., 2012) and in UK and Canada (Clark et al., 2007;Hinton et al., 1997). A likely reason for this is that the DOC dilution primarily occurs during the snowmelt period when large amounts of snowmelt water runoff as overland flow over frozen soil . As we are including events driven by autumn rain and snowmelt, this seasonality difference will not be picked up by the model, but will instead provide a poorer model fit. It can be noted that we have not separated the signal from upland soils from the riparian soils in this study. The riparian zones are included in the signal from both the forest on till soils and forest on sorted sediment. However, till and sorted sediment soils have a different riparian DOC signal. Till soils usually have higher DOC concentrations in the riparian soils (on humid and wet sites) than sediment soil. On humid and wet till soils we have a build-up of riparian peat; hence, water draining to a forest stream on till soils will pass through the riparian peat enriched in DOC and give a higher DOC signal to the forest stream than that of a dry till or a sorted sediment .
In this application of the model, we used a bootstrapping approach to calculate the end-member concentration for the different landscape elements based on headwater stream DOC concentrations. Other approaches for assigning end members could be tested. For example, targeting specific catchments believed to be closer to being true end members can give a better spread of end-member concentrations. Using DOC measurements from soil water from the different landscape elements might also provide more accurate end-member concentration in some circumstances. Grabs et al. (2012) showed that the soil water concentrations from the fine sorted sediments riparian zones are low in Krycklan, averaging around 6 mg L −1 . Furthermore, soil water DOC concentrations in dry till locations have also been found to be relatively low, on average 10 mg L −1 , whereas they were 27 mg L −1 and 33 mg L −1 in humid and wet sites, respectively. By classifying the till-derived soils into three classes (dry, humid, wet) from the topographic wetness index we could potentially improve the predictability of the landscapemixing model for DOC in the catchment. However, this modelling study did not aim to maximise the predictability of DOC in the landscape; instead, the residual analysis can be seen as a diagnostic tool to examine when and where simple land characteristics succeed in explaining the variability in DOC concentration on the landscape scale. Hence, this study should not be seen foremost as a predictive model but rather a learning framework for further development of our conceptual understanding.

Landscape-mixing model performance
The landscape-mixing model (Cooper et al., 2004) combined with the high-resolution DEM offers a simple way of estimating stream DOC concentrations throughout the stream network ( Fig. 1) based on a few headwater observations. Whether the landscape-mixing model is good enough to be used for prediction depends on what the predictions are to be used for, and how much error is acceptable. However, our simple landscape-mixing model produces similar results to the more complicated process-based DOC-3 model (Jutras et al., 2011) where RMSE ranged from 2.4 to 5.1 mg L −1 and R 2 0.27-0.55 (cf. Table 2) in three Nova Scotia streams with stream DOC concentrations 4-40 mg L −1 , similar values to the Krycklan catchment (Fig. 4). Based on the many model measurements calculated (RMSE, RSR, NSE, PBIAS) we assess that the model performed well for two campaigns (May 2003 andSeptember 2008), unsatisfactorily for the one winter campaign (February 2005), and satisfactorily for the remaining four campaigns. How well the model performed depends upon how precisely and accurately the end-member concentrations can be determined (e.g. Fig. 6), whether the modelled solutes behave conservatively, and the degree to which the soil and stream are hydrologically connected (Inamdar et al., 2004). Creed and Band (1998) suggest that stream nutrient export dynamics can be regulated by variable source area dynamics, with large areas of the catchment contributing during high flow but only near-stream zones contributing during periods of low flow. Our study lends support to that idea as the model worked better when the hydrological connectivity to the soils was good, i.e. during high-flow situations (Fig. 4). During baseflow large areas are hydrologically disconnected and hence the landscape-mixing model performed less well (Fig. 4) since it calculates the stream concentrations based on characteristics of the entire catchment. It is also likely that the relatively good model results during high flow are caused by more "conservative-like" behaviour of the DOC due to the shorter in-channel residence times of the stream network (0.5 days) at high flow (Tiwari et al., 2014).

Residual analysis
Using the landscape-mixing model as a baseline, the residual analysis can be used to identify other processes that regulate stream DOC. The residual analysis showed that the model failures were related to hydrological conditions (Fig. 7), indicating that different processes are important for stream DOC during low-vs. high-flow periods. During high/intermediate discharge the landscape-mixing model overestimated DOC in sub-catchments with a high peat coverage (Fig. 8a) while the model underestimated DOC in the same sub-catchments during baseflow (Fig. 8b). Higher concentrations from the wetland during baseflow and lower during high flow would have improved the model according to the residual analysis ( Fig. 8a and b). A possible reason for the inability to predict the peatland-DOC relationship with high accuracy is that the model was constructed on a data set with peat coverage of 0-30 % while it was applied to sub-catchments in the validation data set with a peat coverage of up to 55 % (highlighted as unfilled circles in Fig. 4). Another cause for the model underestimating high values and overestimating low values could be a dampening effect in the bootstrapping approach similar to other regression approaches (Gupta et al., 2009).
The PLS analysis of the residuals showed that during high/intermediate discharge (Fig. 8a) the model underestimates DOC in sub-catchments with much forest and overestimates DOC in sub-catchments with a high peat cover. Given that forests and peatlands are the most common landscape types in the study catchment this makes it difficult to interpret the results, because that means that forest and peatlands are highly correlated (Pearson correlation −0.94; p < 0.001). This can create model artifacts as the overestimated DOC in forest-rich sub-catchments could be because of our inability to capture the true variability in peatland DOC, as discussed above. On the other hand, it could also be a causal relationship. Berggren et al. (2009) showed in a forest-mire gradient in the study area that mixed catchments change their dominant DOC source depending on discharge and that during high discharge the forests become more important as a DOC source.
It should be noted that the landscape-mixing model will only work, in the simple form applied here, if the concentrations downstream are the result of simple mixing of upstream water sources, i.e. if solute transport is conservative. The residual analysis during low-flow situations shows that the mixing model overestimated the DOC concentrations in the lower lying large downstream catchments, with high stream order (Fig. 8b). This highlights the importance of including in-stream processes such as bacterial degradation and/or photo-oxidation of DOC, as well as changing flow paths. These processes need to be included in the conceptual framework when modelling DOC during baseflow. Moody et al. (2013) recorded high photo-oxidation rates (exceeding 10 mg C L −1 day −1 ) during the first 1-2 days of exposure of fresh peat-derived DOC in UK headwater streams, while Köhler et al. (2002) found photo-oxidation rates of the order of 2 mg L −1 day −1 for water from a peat-dominated headwater stream in the Krycklan catchment. In a companion study by Tiwari et al. (2014) we calculated the total photooxidation and bacterial degradation in the Krycklan stream channel network, from headwaters to the outlet, to be less than 1 mg L −1 . Based on this analysis the bacterial degradation and/or photo-oxidation of DOC could only partially explain the overprediction of DOC at downstream sites during low flow.
During baseflow it was the large downstream catchments that had the highest overpredictions of DOC (Fig. 8b). The fact that this landscape type was significant only during baseflow indicates that the overprediction is related to changing flow paths in large catchments during baseflow.  have shown that there is considerable variability in specific discharge in the study catchment and that this affects the DOC exports to the different sites within the catchment. The water in the downstream main stem has a signal more similar to deep groundwater with low DOC concentrations and high base cation concentrations (Klaminder et al., 2011). The overpredicted DOC concentration in the main stem of Krycklan could therefore be related to a large contribution of deeper low DOC groundwater at this scale, diluting the DOC concentrations during baseflow situations. In a companion study by Tiwari et al. (2014) we quantified the amount of deeper groundwater input, using two separate techniques; a mass balance approach and a comparison of specific discharge between a headwater and the outlet. That study showed that at baseflow most of the water (80 %) at the outlet of Krycklan originated from deeper groundwater flow paths, so during low flow the DOC concentration was controlled by the groundwater concentrations and not from the water mixing down from the headwaters. The landscapemixing model assumes that the water at the outlet is the sum of all contributing sources. To create a well-functioning model working during all flow situations one must therefore understand all contributing sources and how they vary during different hydrological situations. By including the deep groundwater as a fourth end member (as a function of discharge) a landscape-mixing model could be constructed that predicts DOC concentrations throughout the stream network during all flow situations (R 2 = 0.88, RMSE = 2.2 mg L −1 ) (Tiwari et al., 2014). This highlights the importance of understanding changing flow paths and seasonal dynamics when modelling DOC in meso-scale catchments.
To conclude, the landscape-mixing model was a useful tool for predicting stream DOC, especially during periods of intermediate and high stream flow. Using the landscapemixing model as a baseline in combination with a residual analysis showed when and where simple mixing did not apply and how the conceptual framework for DOC models must be adapted in space and time. During baseflow it was necessary to consider dilution by low DOC groundwater as another end member in larger downstream catchments.