Cross-scale ensemble projections of dissolved organic carbon dynamics in boreal forest streams

Climate is an important driver of dissolved organic carbon (DOC) dynamics in boreal catchments characterized by networks of streams within forest-wetland landscape mosaics. In this paper, we assess how climate change may affect stream DOC concentrations ([DOC]) and export from boreal forest streams with a multi-model ensemble approach. First, we apply an ensemble of regional climate models (RCMs) to project soil temperatures and stream-flows. These data are then used to drive two biogeochemical models of surface water DOC: (1) The Integrated Catchment model for Carbon (INCA-C), a detailed process-based model of DOC operating at the catchment scale, and (2) The Riparian Integration Model (RIM), a simple dynamic hillslope scale model of stream [DOC]. All RCMs project a consistent increase in temperature and precipitation as well as a shift in spring runoff peaks from May to April. However, they present a considerable range of possible future runoff conditions with an ensemble median increase of 31 % between current and future (2061–2090) conditions. Both biogeochemical models perform well in describing the dynamics of present-day stream [DOC] and fluxes, but disagree in their future projections. Here, we assess possible futures in three boreal catchments representative of forest, mire and mixed landscape elements. INCA-C projects a wider range of stream [DOC] due to its temperature sensitivity, whereas RIM gives consistently larger inter-annual variation and a wider range of exports due to its sensitivity to hydrological variations. The uncertainties associated with modeling complex processes that control future DOC dynamics in boreal and temperate catchments are still the main limitation to our understanding of DOC mechanisms under changing climate conditions. Novel, currently overlooked or unknown drivers may appear that will present new challenges to modelling DOC in the future.


Introduction
Boreal forest landscapes comprise a complex mosaic of forest and peaty wetland elements characterized by organically rich soils that are traversed by a network of streams (Pastor et al. 2003;Laudon et al. 2011). Carbon is mainly exported out of boreal forest catchments in the form of dissolved organic carbon (DOC) (Pastor et al. 2003), which is a key water-quality parameter that influences the fate and dynamics of other carbon driven pollutants in the aquatic environment (Brooks et al. 2007;Kalbitz et al. 2000;Landre et al. 2009). DOC can also act as a medium of carbon exchange in the form of CO 2 emissions from aquatic conduits in boreal systems (Wallin et al. 2013). It is a matter of concern that DOC concentrations ([DOC]) in many parts of the world have been increasing in surface waters Eimers et al. 2008;Oni et al. 2013). Several factors have been identified as drivers of these increases, including declines in acid deposition (Monteith et al. 2007) and landscape disturbances (Schelker et al. 2012) but no scientific consensus has, as yet, been reached because the dominant drivers seem to differ depending on whether local or regional scales of observation are used. For example, recovery from acidification was the main driver of the increase in DOC observed in southern Fennoscandia, the UK and parts of North America Monteith et al. 2007). However, in boreal and other high latitude regions which have not been significantly affected by acidification, there is increasing recognition of the role of climate variability and change as a long-term control of [DOC] (Couture et al. 2012;Laudon et al. 2012Laudon et al. , 2013aLepistö et al. 2014;Tetzlaff et al. 2013;Pastor et al. 2003). Such trends of increasing [DOC] might continue as the future climate warms and precipitation patterns shift in an unpredictable manner (IPCC 2007).
Understanding how climatic factors affect future trajectories of DOC requires an analysis of long-term datasets (Oni et al. 2013), combined with new experimental data (Cleveland et al. 2004) and the development of simulation models (Futter et al. 2007;Wu et al. 2013). While numerous biogeochemical models have been developed to aid our understanding of how climate affects DOC, all contain inherent uncertainties because of our limited understanding of the dominant mechanisms that drive processes at local or regional scales, and the ways in which such processes are represented within models. Our poor understanding of the holistic mechanisms controlling DOC dynamics at a much larger or global scale is also reflected in decisions made about which processes to include or exclude from the model framework.
While no single model can represent all known processes in an all-encompassing way, each can be useful when representing present-day conditions, depending on the scope of the questions the modelers are asking. For example, the simple mixing model of Boyer et al. (2000) focuses mainly on hydrological factors controlling DOC, while the dynamic lacustrine DOC models of Hanson et al. (2004) and Canham et al. (2004) are useful tools for exploring hypotheses concerning DOC dynamics in lakes. The simple model of Jutras et al. (2011) also represents the present-day DOC dynamics well, but is based only on topographic and hydroclimatic factors. Similarly, both DyDOC (Michalzik et al. 2003) and the DOC model proposed by Neff and Asner (2001) are explicit models of processes in the soil phase, which have been useful for predicting soil solution carbon, but their lack of detailed processes occurring in streams might undermine their applicability outside this range. The strength of the model proposed by Lumsdon et al. (2005) is its focus on organic carbon solubility in the soil. The Mire DOC model proposed by Yurova et al. (2008), works well for peaty wetlands, but might not be able to give a credible representation of the contrasting mechanisms controlling DOC in forest streams. The wide variation in process understanding, conceptualization and representation in these DOC models could lead to large uncertainties that could potentially be of the same order of magnitude as the uncertainty associated with future climate projections or even larger. (Larsen et al. 2011). It is therefore important that we understand how these differences can propagate as they may be negligible when simulating present-day conditions but generate significant uncertainties in projections of possible futures.
In the present study we chose two widely differing models of stream DOC. The simple dynamic riparian flow-concentration integration model (RIM) is a hillslope scale model of riparian controls on stream [DOC] in forest catchments ). The concept behind the RIM development was based on the assumption that the riparian zone closest to the stream is a first-order factor controlling stream DOC (Winterdahl et al. 2011a, b) and other streamwater chemistry in forest landscapes (Ledesma et al. 2013). The RIM conceptualization was justified by the fact that the riparian zone is a biogeochemical hotspot where organic-rich soil consistently interacts with ground and surface waters (Grabs et al. 2012) and thus riparian soil water [DOC] and flow paths effectively control stream [DOC]. On the other hand, the Integrated Catchment model for Carbon (INCA-C) is a more complex process-based model of DOC that assumes that both terrestrial and aquatic processes control streamwater DOC at the catchment scale. Combining these models in the form of a process ensemble for the credible assessment of catchment DOC dynamics is therefore a significant challenge in biogeochemistry. Understanding how differences in model assumptions, or the structural representation of processes could amplify the predictive uncertainty of catchment DOC can provide insights that can help in directing us to a time in which arguments over the perceived relative superiority of various current models would be less important, as the future will decide which are best. In this work we conduct a robust assessment of the potential effects of climate change on stream [DOC] and export in boreal forests using an ensemble of regional climate models (RCM) to project soil temperature and stream-flow which, in turn, are used to drive both RIM and INCA-C.

Catchment description
As a case study, we based the cross-scale ensemble modeling of DOC in Svartberget (64°16 0 N, 19°46 0 E), which is a well-monitored and pristine boreal headwater catchment nested at the center of the Krycklan basin in northern Sweden (Laudon et al. 2013b, Fig. 1). However, the approach is not limited specifically to this small catchment since the concept is also applicable at larger scales throughout the boreal or temperate regions. The Svartberget catchment drains an area of about 50 ha and is characterized by forest, peaty wetlands (mire) and mixed forest and mire landscape elements that are typically present in boreal landscapes (Apps et al. 1993;Pastor et al. 2003). The study site is thus a suitable modeling platform on which to investigate the effects of landscape factors on boreal surface water quality and their interactions under future climatic conditions (Fig. 2). Three sites in Svartberget were modelled. The western headwater tributary (denoted as C2 in Fig. 1) drains a completely forested landscape (hereafter refer to as 'forest') while the eastern headwater tributary (denoted as C4 in Fig. 1) drains a mire system (hereafter refer to as 'wetland'). Downstream (denoted as C7 in Fig. 1) is a confluence that drains both streams (hereafter refer to as 'mixed catchment'). The wetland is minerogenic and is covered with Sphagnum spp. (Yurova et al. 2008). The DOC processes in this catchment, therefore, contrast with the forest. Forests in Svartberget consist of century old stands of Norway spruce (Picea abies) and Scots pine (Pinus sylvestris) with some forest understory such as bilberry (Vaccinicium myrtillus) and cowberry (Vaccinicium vitisidaea) (Petrone et al. 2007). The vegetation is in an almost pristine state.
The long-term mean air temperature  in the Svartberget catchment is 1.7°C with minimum and maximum temperatures being -9.5 ± 1.4°C in January and 14.5 ± 1.7°C in July respectively, which provides a longterm mean growing season of *148 days (Oni et al. 2013). Long-term mean annual precipitation is 610 ± 109 mm with snow cover lasting for *170 days (Köhler et al. 2008). Between 35 and 50 % of total annual precipitation falls as snow and winter in the catchment is long (October-May) relative to the short summer period. The mean annual runoff is *320 ± 97 mm as measured using a 90°V-notch weir located at the downstream point C7 (Fig. 1).
The catchment elevation ranges from 235 to 310 m asl (above mean sea level) and is underlain by granite and Gneissic bedrock. Podzol soils are underlain by glacial till in the forest catchment (Berggren et al. 2009) while histosols are common around the wetland catchment ). Svartberget is rich in organic soils that are about 50 cm deep in some parts of the riparian systems in close proximity to the stream (Grabs et al. 2012). Catchment boundaries were delineated using a combination of field observations, light detection and ranging (LIDAR) and a Digital Elevation Model (DEM) (Laudon et al. 2011).

Monitoring data
As the catchment is part of the Svartberget Research Forest, continuous daily temperature and precipitation measurements were available for the period 1981-2010 from a weather station located near the downstream point C7. The downstream specific discharge from C7 was also used for both the forest and wetland catchments. Stream [DOC] was monitored in the forest, wetland and mixed catchments (  intervals. Samples were collected in acid-washed polyethylene bottles and stored frozen before being transported to the laboratory for analysis. Water samples were analyzed in the lab using a Shimadzu 5000 TOC analyzer following common lab protocols as described elsewhere (Köhler et al. 2008;Laudon et al. 2011Winterdahl et al. 2011a;Oni et al. 2013 (Ehret et al. 2012;Teutschbein and Seibert 2013), and do not agree well with observed time series (Fig. 3a, b), these simulations were not used directly as input to the hydrological model in the present study (Fig. 2). A 'distribution mapping' procedure was employed to bias-correct the RCM-simulated precipitation and temperature series on a monthly basis. Distribution mapping, which is also known as 'probability mapping' (Block et al. 2009), 'quantile mapping' (Chen et al. 2013), 'intensity-based statistical downscaling' (Piani et al. 2010), 'histogram equalization' Fig. 2 Schematic representation of the model chain, showing crossscale connectivity and process conceptualization between the Hydrologiska Byråns Vattenbalansavdelning (HBV) rainfall-runoff model , the Integrated Catchment model for Carbon (INCA-C) (Futter et al. 2007), and the Riparian flowconcentration Integrated model (RIM) ). Hydrologically effective rainfall (HER) represents the available water from rainfall and snowmelt that eventually contributes to runoff. Soil moisture deficit (SMD) represents an index of soil dryness derived from the difference between maximum field capacity of a soil and its available soil moisture content. The dotted oval symbolizes the uncertainty envelope surrounding the representation of present-day processes and the subsequent amplification of these uncertainties in potential future conditions moving from hydrologic projections to biogeochemical modelling of streamwater [DOC] and export at a catchment scale (Rojas et al. 2011) or 'distribution-based scaling' (Olsson et al. 2011), has previously been found to be the best correction method for small and meso-scale catchments in Sweden for current (Teutschbein and Seibert 2012) and future climate conditions (Teutschbein and Seibert 2013). The underlying idea is to adjust the theoretical cumulative distribution function (CDF) of RCM-simulated climate values so that it matches the observed CDF.
Several theoretical distributions can be used to describe the probability functions of precipitation and temperature. The Gamma distribution with shape parameter a and scale parameter b is often assumed to be suitable for distributions of precipitation events (Watterson and Dix 2003;Piani et al. 2010), whereas the Gaussian distribution with location parameter l and scale parameter r is usually assumed to fit temperature time-series best (Schoenau and Kehrig 1990). The applied distribution mapping procedure is briefly described below; for a more detailed mathematical description we refer the reader to Teutschbein and Seibert (2012).
As some RCMs tend to simulate a large number of days with low precipitation (i.e., 'drizzle days') instead of being dry, the first step comprised the introduction of an RCMspecific precipitation threshold to avoid substantial distortion of the distribution. In a second step, the distribution parameters were computed for both the observations  and the RCM-simulated control run. They were then used to adjust the RCM-simulated control run climate variables according to Eq. 1, where C is the climate variable of interest (precipitation or temperature); C* represents the bias-corrected climate variable; F stands for the theoretical CDF (Gamma or Gaussian) and F -1 for its inverse; p1 and p2 are the distribution parameters (a and b for Gamma distribution, l and r for Gaussian distribution); the subscripted expression obs indicates observations, and contr stands for the RCM-simulated control run .
In the third and final step, the same distribution parameters (p1 obs , p2 obs , p1 contr and p2 contr ) were applied to adjust the RCM-simulated future scenario run (2061-2090) climate variables (C scen ) according to Eq. 2, where the subscripted expression scen indicates the RCMsimulated scenario run . This final step highlights the underlying assumption of this procedure: the biases are assumed to be stationary and do not changes over time, i.e., the same correction algorithm (distribution parameters) apply to both current and future climate conditions.

Hydrological modeling
The proper representation of catchment hydrology is important for projecting stream [DOC] and fluxes at a catchment scale. The HBV rainfall-runoff model version of Seibert and McDonnell (2010) was employed for hydrological simulations in this study. HBV is a physically based rainfall-runoff model that has been used throughout the world under different hydrological conditions. In the present study, the role of HBV was twofold: first, to represent the hydrological conditions of the catchment under present-day and future conditions; secondly, to provide the necessary hydrological variables (hydrologically effective rainfall and soil moisture deficit) needed by INCA-C (Fig. 2). Running HBV requires daily time-series of air temperature and precipitation. To calibrate the model, temperature, precipitation and runoff for the period 2005-2008 from C7 ( Fig. 1) were used. In the first stage of model calibration, HBV parameters were fine-tuned so that the simulated runoff matched the observed data. This initial parameter set was used as the basis for identifying parameter ranges to use in subsequent Monte Carlo analysis in which the model was run 50,000 times. The parameter set with the highest Nash-Sutcliffe (NS) statistics (Nash and Sutcliffe, 1970) was selected to predict flow from temperature and precipitation. The NS statistic ranges from -? to 1 where a value of 1 indicates perfect agreement between modelled and observed values and a value of 0 indicates that model performance is analogous to replacing modelled values with the mean value of the observed series (Nash and Sutcliffe 1970). Bias-corrected temperature and precipitation series from the 15 RCMs (Table 1) were used as driving input variables to project possible runoff conditions of the catchment for the 2061-2090 time window. The key assumption in HBV projections is that the calibrated parameter values still hold in the future.

Dynamic riparian flow-concentration integration model
The riparian flow-concentration integration model (RIM) is a simple dynamic model that adopts the use of the flowpath controls on streamwater [DOC] along a riparian hillslope gradient of forest catchments (Seibert et al. 2009, Fig. 2). A recent study by Ledesma et al. (2013) has shown that the concept works well for other aspects of stream chemistry such as base cations. The model is driven by flow, stream chemistry, and soil temperature data. The conceptual framework of this approach was therefore based on the flow-related changes on the export of soil water DOC to the stream by the superimposition of lateral fluxes of groundwater through the riparian soil ). The key assumptions of this model are: (1) that all flows go through the organic rich soils in the riparian system, and these flows are therefore the main driver of stream DOC dynamics; (2) that different behaviors of landscape elements within a catchment are reflected within the model without any spatial variability. The static RIM structure originally proposed by Seibert et al. (2009) did not include seasonal influences on DOC dynamics. Consequently, Winterdahl et al. (2011a, b) evaluated different representations of soil temperature factor in the model structure so as to improve seasonal DOC simulations. Winterdahl et al. (2011a) presented a model in which temperature was assumed to have a linear effect on parameter values where the effective parameter value (k 0 ) is a combination of a base value (k 0 ) and a temperature dependent offset (k 1 T): k 0 = k 0 ? k 1 T. In Winterdahl et al. (2011b), soil temperature controls were modelled using an  Table 2.
where T soil denotes soil temperature which was simulated using the model of Rankinen et al. (2004) that is included in the INCA-C terrestrial module (Fig. 2), C base is an empirically estimated base concentration, Q is observed flow, Q 0 is an offset needed to deal with days on which HBV simulated 0 flow and f base represents the base flow:concentration power relationship. RIM was calibrated for the 2006-2010 in each catchment using the ''Solver'' routine in MS-Excel to minimize the sum of squared differences between observed and simulated stream [DOC]. Parameter values from simulations with the highest NS statistics were used for future [DOC] and flux projections driven by bias-corrected climate series.

Process-based integrated model for carbon
The integrated model for carbon (INCA-C) is a processbased biogeochemical model suitable for catchment-scale simulations of the daily dynamics in soil and streamwater [DOC] (Futter et al. 2007(Futter et al. , 2009). The model is best described as a dynamic, semi-distributed model of carbon that integrates terrestrial and aquatic processes to simulate stream [DOC] at a catchment scale (Fig. 2). It has been used at a range of spatial scales including headwater catchments (Futter et al. 2007(Futter et al. , 2009) and larger rivers (Oni et al. 2011;Ledesma et al. 2012) and has been used for simulations based on nearly 50 years of measured data (Lepistö et al. 2014). INCA-C has been widely used for projecting future [DOC] in boreal and temperate surface waters (Aherne et al. 2008;Futter et al. 2009;Oni et al. 2012;Holmberg et al. 2014 inter alia). When using INCA-C, it is necessary to specify several catchment attributes including land cover proportion, spatial boundaries, and length of stream segments, amongst others. The terrestrial interface routes for water travelling through soils and biogeochemical processes were represented for forest and wetland land cover types in this terrestrial module. Further details are provided in Futter et al. (2007).
To run INCA-C, daily time series of air temperature and precipitation from local weather stations were required. These were also used to drive the HBV rainfall-runoff model to generate the necessary hydrological inputs to the INCA-C terrestrial module (Fig. 2). In this case, HBV was used to simulate daily values of hydrologically effective rainfall (HER) and soil moisture deficits (SMD). HER is that fraction of the precipitation that constitutes runoff and contributes to the moisture needed for biogeochemical processes operating within the organic-mineral soil profile continuum (Fig. 2). SMD is an index of soil dryness estimated as the difference between available soil moisture and the maximum water holding capacity of the soil. It is used in INCA-C to set the rate at which moisture-dependent biogeochemical processes in the soil can proceed as described in Futter et al. (2007).
In its terrestrial module, the model also requires estimates of the above-ground litter and soil organic carbon (SOC) pool. Structurally, INCA-C can simulate carbon processes in the soil, both vertically and laterally, by combining the edaphic processes and hydrological flow paths with aquatic processes (Fig. 2). These are represented in a series of interconnected first-order differential equations, which have been fully described in Futter et al. (2007). DOC in solution in the soil can be a product of mineralization of SOC, or derive from root exudates or litter decomposition (Fig. 2). The hydrological connectivity that exists between the terrestrial and aquatic systems is simulated by routing water and carbon from the soil boxes to the stream. The rates of edaphic processes and carbon pools were made dependent on both soil moisture and temperature in order to simulate carbon transformation and production in the terrestrial box. The latter was simulated in INCA-C using a soil temperature model developed by Rankinen et al. (2004). DOC mineralization rates in the stream are dependent on both water temperature and solar radiation. This dependence on temperature and moisture makes the INCA-C model suitable for projecting potential future DOC trajectories at a catchment scale, with the additional possibility of simulating dissolved inorganic carbon (DIC) that can contribute to CO 2 emissions from surface streams (Oni et al. 2012). INCA-C was calibrated using data over the period from 2006 to 2010 in the forest, wetland, and mixed catchments, using a Monte Carlo analysis framework ) in order to more effectively sample the parameter spaces for credible representation of present-day conditions. Model goodness of fit was assessed using NS statistics for simulated flow and in-stream [DOC]. When simulating possible future trajectories of stream DOC, the model was driven by bias-corrected climate series and corresponding hydrological variables.

Precipitation and temperature
Individual uncorrected RCM simulations of temperature (Fig. 3a) and precipitation (Fig. 3b) during the control period 1981-2010 are not in good agreement with observed values. The ensemble median fits the temperature observations better than individual RCMs (Fig. 3a). For precipitation, however, the median is not able to reproduce observed values (Fig. 3b). After applying the distribution mapping procedure to correct for biases, the RCM ensemble median matches the observations during the control period (Fig. 3c, d). Bias-corrected RCM simulations of the future period (2061-2090) project an increase in both temperature (Fig. 3e) and precipitation (Fig. 3f). Overall, temperature is projected to increase between 2.8 and 5.0°C with the median at 3.7°C. Precipitation was projected to increase between 2 and 27 % with the median at 17 %. In terms of seasonal changes, precipitation and temperature are projected to increase for all months of the year (Fig. 3e, f), although these changes are more pronounced during colder months (November-April). The period with temperatures below 0°C was projected to shorten considerably with potential consequences for winter snow accumulation, growing season and spring floods.

Catchment runoff
The HBV rainfall-runoff model performed well in representing the present-day hydrological conditions (NS = 0.82). Although the model slightly underestimated the spring melt peaks, there was good agreement between observed and simulated runoff in winter and summer months. The model overestimated runoff in late autumn/ early winter (Fig. 4a). The multi-RCM projections showed a wide range of future runoff conditions on seasonal (Fig. 4b) and annual scales (Table 3). Most RCMs projected that the spring runoff peak would decrease in the future, while a few models projected slightly higher spring melts than those currently observed. This suggested possible runoff extremes in the future. Despite their differing projections, all 15 RCMs agreed that spring runoff peaks would occur in April, earlier than the present-day peak that occurs in May. Most RCM scenarios projected more runoff in the winter and summer months ( Fig. 4b; Table 3).  Future runoff is projected to be between 305 and 455 mm (median = 392 mm). These values encompass the wide spectrum of projected increases in runoff that range from 2 to 52 % with ensemble median value of 31 % (Table 3).

Present-day DOC simulations
Both INCA-C and RIM performed well in simulating present-day stream [DOC] (Fig. 5). The RIM performed best in the forest catchment (NS = 0.62), and worst in the wetland (NS = 0.52), whereas the mixed catchment resulted in an NS value of 0.54. The RIM performance in the forest catchment is not surprising as the model conceptualization was based on the riparian systems in forest catchments. This means that the RIM performed better in capturing the extreme stream [DOC] values associated with the riparian spring DOC pulses in forest and mixed catchments, but missed some of the low DOC values associated with spring dilution in the wetland. INCA-C also performed best in the forest catchment (NS = 0.52) but its performance range from NS values of 0.49-0.5 in wetland and mixed catchments respectively. This demonstrates the strength of INCA-C in integrating different landscapes, and the edaphic and aquatic processes that are thought to control DOC dynamics at the catchment scale. However, INCA-C missed most of the high spring [DOC] values in the forest and mixed catchments as the model lacked those explicit mechanistic details responsible for driving the riparian spring DOC pulse in this region; these being one of the strengths of the RIM. However, INCA-C performed better in capturing their respective base-flow [DOC] (Fig. 5). This may make INCA-C a better tool for representing how fundamental catchment processes function in reality.

Concentrations
Despite the similar performances of the RIM and INCA-C models in representing present-day stream [DOC], they differed in their projections of future stream [DOC] across the three catchments (Fig. 6). In the forest catchment, all RCMs projected an increase in forest stream [DOC] in RIM, with a median annual mean [DOC] of 19.1 mg/L. INCA-C mostly projected an increased [DOC] with a median annual mean value of 18.1 mg/L, although a decrease was projected with some RCMs. In each case, [DOC] to range between 16.6 and 20.7 mg/L, while RIM projected a narrower range of 17.7-20.9 mg/L ( Fig. 6a; Table 4). However, the inter-annual variation of stream [DOC] generated by RIM appeared to be larger than that generated by INCA-C, despite both models being driven by the same RCM forcings and employing similar assumptions for their calibration strategies.
In the wetland catchment, model outputs exhibited contrasting DOC behaviors from the forest stream (Fig. 6b) In the mixed catchment, both INCA-C and RIM projected similar patterns of stream [DOC] close to the forest catchment, but with a visible mixing effect due to the upland wetland complex (Fig. 6c) (Fig. 6c). This is equivalent to an increase of about 21 % according to the RIM, compared with the decrease (-16 %) and increase (?11 %) projected by INCA-C (Table 4). A consistent pattern of the RIM generating larger inter-annual variation than the INCA-C was also observed in the mixed catchment.

Fluxes
When the projected effects on stream [DOC] from both RIM and INCA-C (Fig. 6) were combined with the runoff projections (Fig. 4) in the form of a DOC flux ensemble (Fig. 2), patterns were obtained that differed from those relating to concentrations (Fig. 7). While both models agreed better on present-day DOC fluxes than on concentrations, they differed significantly with regards to future DOC fluxes and concentrations (Fig. 7). In the forest stream (Fig. 7a), RIM projected a median value for the mean annual DOC flux of 9.1 g/m 2 /year, while INCA-C projected 7.2 g/m 2 /year. In the wetland catchment, the RIM projection for median annual mean DOC flux of 14.9 g/m 2 /year exceeded the 9.0 g/m 2 /year of INCA-C. A similar pattern was obtained in the mixed catchment where the RIM projected a median annual mean DOC flux of 11.2 g/m 2 /year compared with that of 8.8 g/m 2 /year by INCA-C. The ranges of DOC flux projections by all RCMs differed between the RIM and INCA-C models in contrast to their respective projections of concentrations. The RIM projected a wider flux range of 6.4-11.1 g/m 2 /year in the forest catchment, compared with the narrower range of 7.0-8.4 g/m 2 /year projected by INCA-C. In the wetland catchment, the RIM projection of flux ranged from 10.6 to 19.8 g/m 2 /year, while INCA-C projected a range of 8.0-9.9 g/m 2 /year. A similar pattern was obtained in the mixed catchment with RIM projecting a flux range of 8.0-14.0 g/m 2 /year, while INCA-C projected a range of 8.1-9.5 g/m 2 /year. Thus, the RIM projected DOC fluxes to change between ?6 and ?84 % (with the median at 50 %) in forest catchment; by -14 to ?62 % (median = 22 %) in the wetland; and by -2 to ?72 % (median = 38 %) in the mixed catchment ( Fig. 7; Table 4). This is in contrast to the narrower ranges projected by INCA-C; ?12 to ?34 % in forest, -34 to -19 % in the wetland, and -4 to ?13 % in the mixed catchment with ensemble median of 24, -26 and ?5 % respectively. RIM projections consistently showed larger inter-annual variations in all catchments. The inter-annual variation in DOC fluxes projected by each RCM in INCA-C was lower in wetland than in either the forest or mixed catchments.

Change in hydroclimatic regimes
A growing number of studies show that climate change will greatly affect hydroclimatic regimes and DOC exports in environments at high latitudes if future conditions are warmer and wetter (Jungqvist et al. 2014;Carey et al. 2010;Tetzlaff et al. 2013;Mellander et al. 2007;Lepistö et al. 2014). This study supports these earlier findings, with the RCM ensemble projecting a wide range of possible future hydroclimatic conditions (Table 3). Precipitation is projected to increase 17 % and temperature 3.7°C (based on the ensemble median). At annual (Table 3) and seasonal ( Fig. 4) scales, runoff conditions were projected to cover a wider range than that projected for precipitation (Table 3), and the ensemble mean suggests that spring peaks could decline by up to 13 % by the end of the present century (Fig. 4b).
Model projections indicated that extreme values could be more pronounced in the future. For example, the annual runoff of individual years ranges between 128 and 576 mm for current conditions, but can range between 123 and 867 mm for future conditions. A recent study by Chou et al. (2013) demonstrated a similar possibility that the range of precipitation extremes might increase in the future to such an extent that they could intensify dry and wet episodes. They also argued that the seasonal cycle would have pronounced effects on hydroclimatic extremes, even if annual rainfall were to be relatively stable in the future.
Any uncertainties in current HBV simulations of rainfall-runoff conditions could be constrained in future studies by using an ensemble of hydrological models. However, in the present study, all RCMs agreed well on the probability that peak spring runoff events would shift from May to April (Fig. 4b). Such a shift in spring runoff events is an indication of the possible extension of the growing season in the boreal region. Longer growing seasons and warmer temperatures (Jungqvist et al. 2014) could lead to ecological succession, changes in litter production (timing, and quality), and changes in stream chemistry (Euskirchen et al. 2006;Linderholm, 2006;Norby and Cotrufo 1998). This also implies the possibility of a regime shift where winter precipitation patterns are partitioned more toward rainfall dominance (Dore 2005), indicating a more biogeochemically active winter season (Laudon et al. 2013a). While the actual magnitude of such changes in the future hydroclimatic conditions might be missed, the direction of As highlighted above, a change in hydroclimatic regimes could drive biogeochemical processes beyond the range of our current experience and understanding. This would have implications for predictive modelling. Observing differences that occur at different scales and between different processes has added to our knowledge concerning those factors that control DOC dynamics (Clark et al. 2010). The development of models based on a mechanistic understanding has also helped to unify our understanding of processes (Futter et al. 2007;Rowe et al. 2014). In the present study we have shown that even apparently negligible differences in simulating present day conditions can lead to large uncertainties when projecting future catchment [DOC] and exports (Fig. 8). This is because our current model development relies heavily on the perceptual understanding of local processes, the types of data available to model developers, and key processes that are included. All these factors emphasize how the strength of a model can make it a useful tool for addressing pertinent questions, but also how decisions by the model developer about which processes to simulate or exclude may amplify into a large degree of prediction uncertainty (Fig. 8). This is because all models assume that present structures also hold under future conditions; but in reality, climate change could drive biogeochemical processes to such novel states that our present assumptions would be undermined. For example, a 'known unknown' is that a change in hydroclimatic regimes projected by all RCMs could make the dormant seasons more active as the ambient environment favors more biotic or bacterial activity, or as water and DOC fluxes become synchronous (Laudon et al. 2013a). The intensification of biogeochemical processes during any future extension of the growing season could also present other 'unknown unknowns' in the form of novel mechanisms that might drive future DOC dynamics beyond the boundaries of our current level of understanding. The present study was designed as a pilot study with the principal aim of promoting cross-scale approaches to projecting catchment DOC and indicating how future models should be tested and validated. For simplicity, we limited the study to two models, choosing two significantly different models in order to assess where the dominant b Fig. 7 Cross-scale ensemble projections of annual stream DOC fluxes in INCA-C (y-axis) and RIM (x-axis) relative to the presentday conditions in forest (a), wetland (b) and the mixed catchment (c). Standard errors are shown in red (RIM) and dark blue (INCA-C). Symbols represent the mean fluxes of each RCM forcing and the respective RCMs are denoted as in-symbol numbering following their sequence in Table 1. RIM projected median annual mean DOC fluxes of 9.1 versus 7.7 g/m 2 /year in INCA-C (forest), 14.9 versus 9.0 g/m 2 / year (wetland) and 11.2 versus 8.8 g/m 2 /year (mixed) sources of uncertainty are, and how they might be further amplified under climate change. The RIM is a simple dynamic, hillslope focused model, which operates with the main assumption that the near-stream riparian zones in forest catchments are the most important drivers of stream [DOC] and exports. It excludes any consideration of spatial variability that might influence the overall catchment behavior Winterdahl et al. 2011a, b). By contrast, INCA-C is a process-based and semi-distributed model that uses a catchment scale approach to represent stream DOC (Futter et al. 2007;Oni et al. 2012). The key assumption in INCA-C was that whole catchment processes (both in the terrestrial and aquatic compartments) are equally important and should be integrated for credible predictions to be made of stream [DOC] and exports. The detailed process representation in INCA-C led to there being a more parameters to be estimated (Futter et al. 2007) and so could lead to equifinality (Beven 2006). In the work presented here, we have focused mainly on the consequences of biogeochemical model structural uncertainty on possible future [DOC] under an ensemble of RCM climate projections. We have not explored all aspects of the uncertainties associated with equifinality of multiple parameter sets. We have based our hydrological results on a single parameterization of the HBV model. Oni et al. (2014) have noted that more credible projections of future hydrological conditions can be obtained with an ensemble of behavioral parameter sets than with a single realization.
We have also used on a single parameterization of the RIM model at each of the three modelled catchments but it should be noted that Winterdahl et al. (2011a, b) have used an ensemble of RIM parameter sets to describe uncertainty in present day DOC predictions. While a Monte Carlo analysis was conducted to identify an ensemble of plausible INCA-C parameter sets, only a single parameter set was used for the climate change simulations presented here. Using ensembles of parameter sets could help to better constrain the uncertainty surrounding projections of future [DOC] but would be computationally challenging.
The performance of INCA-C is encouraging for a multiparameterized biogeochemical model that requires the simultaneous calibration of several parameters in the terrestrial and aquatic boxes (Fig. 2). However, no single model can be all encompassing since they all have inherent strengths and weaknesses when representing present-day biogeochemical conditions. For example, we observed that INCA-C missed some of the DOC peaks despite its robust representation of landscape controls and its ability to give reasonable simulations of stream DOC when using a credible parameter optimization strategy. This shortcoming can be attributed to a riparian DOC surge that resulted from the continuous inundation of the riparian system as groundwater continually rises to activate the carbon pool in the topsoil layers. Where the INCA-C model has been previously applied, it has performed well in capturing the whole inter-annual range including the peak values (Futter et al. 2007;Oni et al. 2011). This cross-scale approach has helped us to recognize a need for the INCA-C structure to be revised to include a mechanism for generating the riparian DOC pulse that is simulated by RIM under highflow conditions. This modification will make INCA-C a better tool for understanding DOC dynamics and for a robust uncertainty assessment of DOC dynamics under hydroclimatic extremes.
By contrast, the fact that the RIM concept was based on the riparian zone mechanism can explain why this model performed better than INCA-C in capturing the present-day peak spring DOC pulses. Despite its strength in the forest catchment, RIM missed some of the low [DOC] associated with spring dilution in the wetland. This is because of the presence of minerogenic mires with contrasting dilution mechanisms, and upland forests having possible roles that were not duly considered in the current model structure Winterdahl et al. 2011a). For example, Yurova et al. (2008) demonstrated the importance of sorption in modeling DOC outputs from boreal mires, but this was lacking in RIM. As a result RIM missed some of the low spring [DOC] and overestimated the high summer [DOC] values in the wetland. This can make the model structure depart from reality toward the downstream reaches of large catchment as the spatial heterogeneity of landscape control on groundwater dynamics can also make the riparian influence on stream [DOC] vary along a river gradient of meso-scale to large catchments (Tiwari et al. 2014). However, the model structure assumed a constant flow-concentration relationship. This reinforces the potential shortcomings of the RIM model when applied outside of its original conceptual landscape framework. The RCMs are fairly consistent in their projections of a warmer and wetter future but present a considerable range of possible future scenarios. This has implications on how the RIM and INCA-C models propagated stream [DOC] and fluxes into the future. For example, INCA-C consistently projected wider ranges of stream [DOC] but the inter-annual variations projected by RIM were larger. The wider range of stream [DOC] projected by INCA-C might be attributed to its climate forcing being very temperature sensitive. This might be explained by studies that have shown the terrestrial carbon pool to be very temperature sensitive in respect of climate change effects (Davidson and Janssens 2006;Lindroth et al. 1998). This makes INCA-C respond more rapidly to the warmer world projected by all RCMs and might explain the wide range of [DOC] projections. However, INCA-C missed the hydrological sensitivity necessary to capture the riparian DOC pulses, and this might explain why its projected fluxes were considerably lower than those of RIM. This may also be explained by some of the DOC being lost in the terrestrial and aquatic carbon pool (Fig. 2) as a form of positive feedback in response to warmer temperatures (Davidson and Janssens 2006).
The large inter-annual variations observed in the RIM projections may be explained by sensitivity to the hydrological variability that would become more common in a wetter future and due to the lack of any feedback between climate-related controls on DOC production and consumption. The RIM formulation presented here simulates a monotonic relationship between soil temperature and [DOC] whereas INCA-C can simulate a non-monotonic relationship between soil temperature and stream [DOC] as it explicitly includes production and consumption terms. RIM does not have any way of representing long-term changes that could occur due to accumulation or depletion of soil carbon pools, which would affect DOC fluxes. The smaller range in INCA-C-projected [DOC] suggests that a balancing feedback mechanism (such as microbial/photo mineralization, sorption-desorption, or terrestrial soil carbon depletion processes) would become more important in a warmer, wetter future and change the underlying flow-concentration relationship conceptualized in RIM. We also suggest that while apparently negligible differences between model assumptions may lead to similar predictions for present-day conditions, the uncertainties will amplify under future conditions.

Conclusion
We have shown that output from both DOC models agreed well in simulating present-day [DOC] and fluxes despite differences in their structural design and assumptions. This implies that both models are equally useful tools when seeking answers to relevant questions framed under current conditions and that surface water [DOC] can be understood as the product of either catchment or hillslope scale processes. However, future projections from the two models differed. Small or negligible present day differences in model assumptions were amplified by the range of possible future climate scenarios. Combining DOC models in the form of a process ensemble is more difficult than creating climate model projections per se due to the large number of differences in the way different models represent processes or suites of component processes. However, the use of an ensemble of climate models, together with cross-scale comparisons of biogeochemical models applied to a range of landscape element types is an important technique when assessing possible future trajectories of DOC in boreal and temperate streams. For example, INCA-C could not capture the high observed spring DOC pulses resulting from a strong riparian system influence while RIM could simulate spring peaks that INCA-C missed. The significant differences between the two model structures amplify differences in future hydroclimatic sensitivity. RIM projected wider inter-annual variations and range of fluxes due to a complete lack of any balancing feedback mechanism of climate change from other upland processes. Moreover, the amalgamation of whole catchment behavior within the riparian system challenged some of the fundamental assumptions of RIM under future conditions. If the explicit riparian DOC pool represented in the hillslope-scale RIM were included in INCA-C, it might improve its usefulness in the northern boreal region.