A coupled hydrology–biogeochemistry model to simulate dissolved organic carbon exports from a permafrost‐influenced catchment

We outline the development of a simple, coupled hydrology–biogeochemistry model for simulating stream discharge and dissolved organic carbon (DOC) dynamics in data sparse, permafrost‐influenced catchments with large stores of soil organic carbon. The model incorporates the influence of active layer dynamics and slope aspect on hydrological flowpaths and resulting DOC mobilization. Calibration and evaluation of the model was undertaken using observations from Granger Basin within the Wolf Creek research basin, Yukon, northern Canada. Results show that the model was able to capture the dominant hydrological response and DOC dynamics of the catchment reasonably well. Simulated DOC was highly correlated with observed DOC (r2 = 0.65) for the study period. During the snowmelt period, the model adequately captured the observed dynamics, with simulations generally reflecting the timing and magnitude of the observed DOC and stream discharge. The model was less successful over the later summer period although this partly reflected a lack of DOC observations for calibration. The developed model offers a valuable framework for investigating the interactions between hydrological and DOC processes in these highly dynamic systems, where data acquisition is often very difficult. © 2015 The Authors Hydrological Processes Published by John Wiley & Sons, Ltd.


INTRODUCTION
The Arctic and sub-arctic regions are dominated by permafrost and discontinuous permafrost soils, which contain large stores of soil organic carbon (SOC) (Guo et al., 2007;Schuur et al., 2008). Proportionally large amounts of carbon are released from these catchments via aquatic pathways (Cole et al., 2007;Schuur et al., 2008;Kicklighter et al., 2013) with dissolved organic carbon (DOC) being an important form of carbon in these streams (Cole et al., 2007;Balcarczyk et al., 2009;Cory et al., 2014). Permafrost has a significant influence on hydrological flowpaths (Quinton and Marsh, 1999;Gray et al., 2001;Carey and Woo, 2001) with the thawed active layer controlling the timing of DOC production and exports (Carey, 2003;Guo et al., 2012;Kicklighter et al., 2013). Recent work indicates a deepening of the active layer (Oelke et al., 2004;Quinton et al., 2011;Kicklighter et al., 2013) and the degradation and loss of discontinuous permafrost soils from surface temperature warming, which will have large impacts on both the hydrological and DOC cycles (Dye, 2003;Oelke et al., 2004;Finlay et al., 2006;Cole et al., 2007;McGuire and Anderson, 2009;Quinton et al., 2011;Cory et al., 2014).
Studies investigating the interaction between SOC and aquatic pathways in high latitude areas have often focussed on large pan-arctic catchments (Striegl et al., 2005;Guo et al., 2012;Holmes et al., 2012;Mann et al., 2012;Tank et al., 2012). Understandably, these larger catchments have received greater attention as collectively they represent roughly 10% of the earth's discharge and DOC exports (Guo and Macdonald, 2006;Wickland et al., 2012). These large-scale studies have revealed valuable information on the processes regulating DOC exports in these systems and the underlying controlling factors. Studies focusing on controls of DOC at smaller scales have shown that the export of DOC to streams is influenced by the presence and location of permafrost active layer thaw depth and runoff processes within the catchment (Carey, 2003;Kawahigashi et al., 2004;O'Donnell and Jones, 2006;Petrone et al., 2006;Prokushkin et al., 2007;Koch et al., 2013). Runoff generation in high-latitude catchments is strongly influenced by snowmelt processes (Mcnamara et al., 1997;Dornes et al., 2008;Carey and Woo, 2001) with much of the snowmelt entering the stream as overland flow or through the upper organic layer, which is characterized by relatively high infiltration rates (Quinton and Marsh, 1999;Carey and Woo, 2001;Gray et al., 2001;Carey and Quinton, 2005;Koch et al., 2013). Permafrost soil in these regions limits deep drainage, sustaining near surface water tables influencing runoff generation (Mcnamara et al., 1997;Quinton and Marsh, 1999;Carey and Woo, 1999). In addition, the presence of organic soil with high infiltration rates and inter-hummock hollows in the organic layer support dominant near surface flowpaths (Mcnamara et al., 1997;Quinton and Marsh, 1999).
Dissolved organic carbon exports from permafrostdominated catchments are high during early spring when the thawed active layer is shallow (Carey, 2003;Frey and McClelland, 2009) and soil water is in contact with relatively modern labile forms of SOC (Guo and Macdonald, 2006;Raymond et al., 2007). Later in the season, as the active layer deepens, soil water comes into contact with mineral soils associated with less labile carbon and acts to immobilize DOC via sorption (Lyon et al., 2010;Koch et al., 2013). However, observations of older carbon within stream DOC during late summer Guo and Macdonald, 2006) highlight the importance of the contribution of DOC from deeper SOC within mineral soils and often evolved from instream processing (Zou et al., 2006). With projections of a loss of areal coverage of permafrost and a deepening of the active layer in sub-arctic regions (Cole et al., 2007;McClelland et al., 2007), it is crucial to understand the interaction between hydrological flowpaths and DOC production and transport from soils to streams, which ultimately enter and influence biogeochemical processes in the Arctic ocean McGuire and Anderson, 2009;Holmes et al., 2012;Cory et al., 2014).
Several studies have modelled the export of DOC in streams across different climatic regions using a variety of model complexity and structures. Linear models are commonly used to estimate DOC in sub-arctic environments based on a proportional relationship between stream discharge and DOC (Carey, 2003;Jutras et al., 2011;Guo et al., 2012). Using this method of statistically based estimation provides an overall understanding of seasonal or annual exports, but only limited understanding of the underlying controlling processes. An alternative approach is the use of physically based models. Several physically based models have been developed to simulate DOC production in soils and are used in conjunction with hydrological models to estimate DOC catchment exports. This approach combines the output from hydrological models with that of terrestrial based soil carbon models, with the hydrological and biogeochemical models calibrated individually. One of the most widely applied of these is the Dynamic DOC model (Michalzik et al., 2003), which combines soil carbon production and loss functions for multiple soil layers and includes a simple hydrological model to simulate soil moisture and runoff processes. Recently, Jutras et al. (2011) used the Forest Hydrological model to supply soil temperature and moisture to the biogeochemical based DOC-3 model to simulate DOC exports in 11 catchments in Nova Scotia, Canada. Other DOC models include the Integrated Catchments Model for Carbon (Futter et al., 2007) and the Estimating Carbon in Organic Soils Sequestration and Emissions (ECOSSE) (Smith et al., 2010a(Smith et al., ,2010bAbdalla et al., 2014). However, these models aim to represent all organic carbon transformations and, therefore, may be overly complex when trying to simulate DOC dynamics at the headwater scale in high-latitude environments, where data acquisition is often difficult. Alternatively, Birkel et al. (2014) and Dick et al. (2014) present a parsimonious hydrological-biogeochemical model, which incorporates the main controlling processes to simulate DOC and stream discharge on a daily basis from an upland Scottish catchment. Such use of a parsimonious coupled model provides a framework to identify the dominant processes affecting DOC transport at a catchment scale.
Here, we outline the development and evaluation of a parsimonious hydrological-biogeochemical model based on the dominant processes of DOC production and export to simulate daily DOC and stream discharge in a small, but complex, sub-arctic alpine catchment with permafrost influence.
The specific objectives are as follows: • Develop and calibrate a simple, coupled hydrological and biogeochemical model to simulate stream discharge and DOC dynamics in a sub-arctic alpine catchment. • Evaluate the suitability of such coupled models for representing dominant hydrological and biogeochemical processes in these systems.

STUDY SITE AND DATA
The 7.4-km 2 Granger basin is a sub-arctic alpine catchment within the Wolf Creek research basin, Yukon, Canada ( Figure 1), with an elevation range of 1310-2250 m.a.s.l. . Annual precipitation in the catchment is estimated at 478 mm (Laudon et al., 2012) with strong seasonality and~50% falling as snow (Carey, 2003). Air temperature varies with a winter and 5384 J. S. LESSELS ET AL.
summer mean of À21 and 15°C, respectively (Carey, 2003). The soils are characterized by 0.05-m to 0.4-m thick organic layers overlying fine to coarse textured colluvium (Carey and Woo, 1999). Permafrost is dominant on north-facing slopes at higher elevations (Carey, 2003;Lewkowicz and Ednie, 2004) with seasonal frost on south-facing slopes. The snow depth, organic layer thickness and presence of permafrost are all heavily dependent on the slope aspect, with thicker organic layers, and higher likelihood of permafrost and late lying deeper snow packs on north-facing slopes throughout the catchment (Carey, 2003;Lewkowicz and Ednie, 2004;McCartney et al., 2006). An estimated 43% of the greater Wolf Creek research basin is underlain by permafrost (Lewkowicz and Ednie, 2004). Field data used here were collected over a 7-year period with DOC samples collected over 5 years during and after the spring snowmelt period of each year. Stream discharge was measured at the outlet during summer months, with gauging equipment removed during winter months when the stream freezes. A total of 179 DOC stream water samples were collected during the study period. In addition to stream water samples, shallow groundwater DOC samples were collected at several locations along a transect across north-facing and south-facing slopes perpendicular to the main stream with the approximate location of the transect (Figure 3). The shallow groundwater dip-wells were constructed using 35-mm PVC pipe with screen mesh along the entire below surface section to a depth of 0.75 and 1.0 m on the north-facing and south-facing slopes, respectively (Carey, 2003). Sampling of the wells commenced on 9 May 2002 but was later for some wells because of snow cover or lack of water. South-facing slopes averaged 25.8 mg C L À1 between 10 and 15 May and later reduced to an average of 7.8 mg C L À1 between 1 and 5 June while the north-facing slopes averaged 55.4 mg C L À1 between 15 and 20 May and 19.1 mg C L À1 during 1 and 5 June (Carey, 2003). Additional information about the DOC samples is provided by Carey (2003). Air temperature and daily rainfall data were obtained from a nearby weather station (c.a. 2.5 km). To fill gaps of missing precipitation and air temperature, additional data were obtained from the meteorological station at Whitehorse airport, located 15 km from the catchment and at an elevation of 706 m.a.s.l. (Carey, 2003). Missing air temperature data were estimated using a linear relationship between the alpine station and Whitehorse airport, while no correction was undertaken for precipitation data. Daily potential evapotranspiration was estimated using the Thornthwaite equation (Thornthwaite, 1948) corrected for daily temperature using the SPEI package (Beguería and Vicente-Serrano, 2013).

THE COUPLED HYDROLOGY-DOC MODEL
A semi-distributed conceptual model was developed to simulate stream discharge and DOC production and loss. The coupled model consists of two parts: a hydrological component using a modified version of the Hydrologiska Byråns Vattenbalansavdelning (HBV) rainfall-runoff model (Lindström et al., 1997) and a DOC production component using a modified version of the DOC production and loss functions based on the algorithms in the ECOSSE soil carbon model (Smith et al., 2007). Figure 2 presents the conceptualization of the two components and how they are linked together. Both model components have been programmed in the R statistical language (R Core Team, 2012) and C++ using the Rcpp (Eddelbuettel and François, 2011) and the RcppArmadillo packages (Eddelbuettel and Sanderson, 2014).

Hydrological model component
The HBV model was further modified to account for the influence of permafrost and seasonal frost on stream discharge. The HBV model has several subroutines that account for snow accumulation and melt, evaporation, soil moisture and runoff described in detail by (Lindström et al., 1997). A description of each coefficient used in the model is outlined in Table I.
Here, we accounted for the effects of aspect on snow accumulation and the onset of snowmelt. On the basis of the work of Hottelet and Braun (1993), the model was applied as a semi-distributed model, separated into three hydrological response units (HRUs) based on aspect: north, south and east-west units. To minimize the number of model coefficients, the same model coefficients were used for all HRUs but were altered by the aspect coefficient of the model. This method modifies the degree-day factor on south-facing slopes by multiplying the degree-day constant (CFMAX) by FASPECT on south-facing slopes and dividing by FASPECT on north-facing slopes (Hottelet and Braun, 1993;Konz and Seibert, 2010). The overall effect of this modification is that north-facing slopes generally have a larger snowpack and the onset of the snowmelt is later than that of south-facing slopes.
In addition, all water storage components have been modified using a degree-day function to incorporate the effect of permafrost and seasonal frost on the hydrology of the catchment. The standard HBV model has three storage components: soil moisture, and upper and lower groundwater storages. On the basis of the work of Heerma (2013), the same degree-day factor used for snowmelt was used to represent the freeze-thaw processes of the permafrost active layer and seasonal frost in each storage component for all three HRUs. Conceptually, the water stored in the ground gradually freezes when the temperature drops below the temperature threshold of the model based on the following equations: where Y (t) represents the frozen water content of the soil, or upper or lower storage zones. When air temperature is above the threshold, thawing of the three storages commences based on the following equations:  This conceptual thaw-freeze process restricts runoff from the upper and lower storages during winter months. To minimize model complexity, the freezing and thawing of these storages are undertaken independent of each other. The freezing of the soil moisture within the model increases the concentration of stored DOC and also limits the production of DOC from the reduction in soil moisture.
Three linear response functions were used: two runoff components for the upper storage, where the first reflects overland flow when the upper storage is at or greater than the field capacity. The second response function reflects the subsurface drainage through the organic layer. The third response function reflects base-flow runoff from deeper storages.

DOC model component
The DOC component of the coupled model consists of a production and loss routine based on the methods used in the ECOSSE soil carbon model (Smith et al., 2007). The coefficients of the biogeochemical component of the model are detailed in Table II. The DOC production function uses a single carbon pool, does not consider temporal inputs from plant material and assumes a static total organic carbon (TOC) pool. Therefore, TOC is the only carbon pool used in this study, and the production rate constant of Smith et al. (2007) for the different carbon pools was not used but estimated during model calibration. The initial range of the production coefficient (k prod ) was based on the range of the coefficients used for the carbon pools in the ECOSSE model (Table IV). The initial range of the coefficient representing the decomposition to biomass and CO 2 is also provided in Table II. TOC was set at 5.4 kg m 2 for the catchment on the basis of the estimates by Hugelius et al. (2013) for the upper 30 cm of the soil profile representing an average for the catchment.
Daily DOC production was estimated on the basis of a combination of the DOC production coefficient and two additional terms using soil moisture, air temperature and a rate constant modifier, where k prod represents the production rate constant (estimated during calibration), TOC is the total organic carbon (kg m À2 ) in the upper 30 cm of the soil, and the rate modifiers, a and b, are reflective of the air temperature and soil moisture, respectively. Therefore, DOC production was assumed to be dominant in the upper 30 cm of the soil profile similar to Michalzik et al. (2003) and Smith et al. (2007). The temperature rate modifier was a based on a Q 10 temperature coefficient relationship with air temperature where and T (t) is the air temperature at time t. A fixed Q 10 rate of 2 was used to minimize model complexity, which is commonly used in soil carbon models (Tjoelker et al., 2001) and falls within the range of Q 10 rates observed in Arctic soils (Lenton and Huntingford, 2003). With a value of 2, the effect of the temperature rate modifier a (t) doubles for every 10-°C temperature increase. This was also used by Michalzik et al. (2003) for DOC production simulations within the Dynamic DOC model. The rate modifier for soil moisture was based on that of ECOSSE (Smith et al., 2007), but altered for use with the hydrological model rather than the physically based where SM (t) is the soil moisture at time t and FC is the estimated field capacity. The loss of DOC from decomposition into biomass was based on the functions of ECOSSE where and any decomposition of DOC is assumed to be lost to CO 2 or sorbed to the mineral soil.

Model calibration
Model calibration was performed using the Generalized Likelihood Uncertainty Estimate methodology (Beven and Binley, 1992). Many goodness-of-fit measures have limitations; for example, the Nash-Sutcliffe Efficiency measure is the most commonly used likelihood function applied in conjunction with Generalized Likelihood Uncertainty Estimate (Schoups and Vrugt, 2010) but is explicitly biassed towards high flows (Criss and Winston, 2008;Gupta et al., 2009). This is problematic in calibrating hydrological models in alpine catchments with large snowmelt flows and long periods of winter base-flow conditions (Hamilton et al., 2000). Therefore,  (Criss and Winston, 2008) was used instead of the Nash-Sutcliffe Efficiency to minimize the effect of large flows on model performance. In addition, the coefficient of determination (r 2 ) multiplied by the slope of the regression (Br 2 ) (Krause et al., 2005) between the observed and the predicted values was also used to evaluate the model efficiency for stream discharge. The use of Br 2 was chosen in this study as it was found to provide a more robust representation of the relationship between the simulated and observed values and less susceptible to larger observed values. DOC simulations were also evaluated by Br 2 as this provided a more robust measure compared with the r 2 . The three model efficiency measures were combined to formulate a single objective function, We based calibration on hydrological data from 2001 to 2003 and DOC observations collected during the summers of 2002 and 2003. Because of a limited amount of available data during winter flows when the stream freezes, calibration of stream discharge and DOC was restricted to summer periods where daily stream discharge observations were made. As the TOC and hydrological storages of the model are dynamic in nature, a spin-up period of 40 years was used to initialize the model. This long spin-up period was chosen to ensure the model was in a steady state prior to the study period. A Monte Carlo approach was then used to evaluate a suitable range of model coefficients using 10 6 model evaluations. The behavioural parameter set chosen on the basis of the 0.95 percentile of the objective function and corresponded to the best 50 000 parameter sets.
To investigate the sensitivity and influence of the model coefficients on the simulations, a global sensitivity analysis was conducted. The sensitivity analysis was undertaken using the Latin Hypercube one-factor-at-atime method proposed by van Griensven et al. (2006) as programmed in the hydroPSO R package (Zambrano-Bigiarini and Rojas, 2013). This procedure uses an efficient Latin Hypercube sampling technique to evaluate the model efficiency over the range of the model coefficients to provide a rank of importance for each model coefficient. Figure 3 shows the observed stream discharge and DOC concentrations for the period between 2001 and 2008. Stream discharge is highest during the snowmelt period and then decreases throughout summer; this trend continues until early winter. During this period, stream discharge is responsive to large rainfall events, which results in large discharge peaks later in the season (e.g. 2006 and 2007). Table III provides a summary of annual statistics of stream discharge. Minimum daily discharge is less than 1 mm day À1 for each year, and mean discharge is less than 2 mm day À1 . There are large inter-annual variations in maximum stream discharge with a range of 5.2-25.7 mm day À1 across the study period. As with  discharge, the highest DOC concentrations are closely linked to the snowmelt runoff period with the observed DOC peak occurring during the snowmelt event before but prior to the peak of stream discharge. DOC concentrations vary from 0.5 to 25.5 mg L À1 throughout the study period (Table III). Of the 6 years of DOC observations, only one year (2003) did not have any samples greater than 20 mg L À1 . Table V provides the initial range and estimated posterior mean and range for each model coefficient. In addition, the table also provides an estimate of the sensitivity of each coefficient. Of the hydrological model coefficients, the snowfall correction factor (SFCF) was found to be the most influential to the model simulations. Hydrological runoff coefficients k1 and UZL were the least important coefficients in the model, indicating that the quick runoff component of the model has little influence on the hydrological responses in the model. Figure 4 shows the simulated stream discharge for the calibration period (2001)(2002)(2003) with 90% uncertainty bands. The simulations capture the observed stream discharge dynamics in each year quite well although there are differences. Simulation results are satisfactory, considering the simplicity (low parameter numbers) of the model applied, with no systematic overestimations or underestimations. The simulated discharge during 2001 matches the observed onset of the spring melt event adequately, but the mean simulated maximum peak is not as high as the observed peak. Observed peak discharge during the spring melt is 26 mm day À1 and does fall  Figure 5 presents simulated and observed stream discharges during the years 2006-2008 to test the model on data not used in the calibration. The simulated discharge dynamics closely follow that of the observed stream discharge for each year. Onset and magnitude of the spring melt event of the simulated discharge closely resembles that of the observed discharge for all 3 years. During this period, the simulated discharge peaks resulting from rainfall events later in the season closely resemble the observed stream discharge, with all discharge observations falling within the uncertainty bounds of the simulations.

Simulation of DOC dynamics
For the calibration period, in 2002, the mean simulated DOC shows two concentration peaks coinciding with the spring melt event (Figure 6). There are no observed DOC samples during the first simulated DOC peak, and on the basis of direct field observations, the first sampled DOC during this year captured the spring melt event. The simulated DOC is similar to that of the observed DOC, but occurs before the observed DOC event, possibly due to the small simulated stream discharge event prior to the actual observed snowmelt event, which was most likely caused by a short period of above-freezing air temperatures. In contrast, the simulated DOC for 2003 closely resembles the dynamics of the observed with the majority of observations falling within the bounds of the simulated uncertainty. Figure 7 shows the model tests for the years not used in calibration. Simulations generally follow the dynamics of observed DOC for all 3 years. The commencement of the simulated exports of the spring event matches that of the observed DOC. The duration of the simulated spring event is also similar to the observed, but the maximum simulated DOC is consistently below that observed. However, the majority of observed DOC samples fall within the simulated uncertainty bounds for all 3 years. For all 5 years (including the calibration period), simulated DOC concentrations show high variability associated with smaller runoff events during the end of the summer period, which were not captured during field campaigns. Uncertainty bounds are greatest during the recession of the spring melt event and later in the season where there are very few observations. Dissolved organic carbon production in the model occurs within the soil storage ( Figure 2) and is controlled by soil water content and air temperature. Figure 8 shows the mean simulated DOC for the south-facing and northfacing slopes with observed DOC samples collected from dip wells during 2002. Simulated DOC accumulates during the winter period as the amount of frozen soil water increases and is higher at the south-facing slope. The rapid reduction in the concentration of soil water DOC in both slopes shown in the observed dipwell samples is well captured with the simulations. After this rapid flushing of soil water DOC, the concentration of DOC in the soil water increases again over the remaining summer months. During the snowmelt event, DOC is observed to reach a maximum concentration during the rising limb of the freshet hydrograph. Figure 9 shows the simulated and observed stream runoff and DOC exports during 2008. The simulations of the model closely resemble the dynamics of the observed stream discharge and DOC. Simulated DOC peaked during the snowmelt period on the rising limb of the stream discharge.
Overall, the model simulated DOC relatively well during the study period, but there was a wide range of model performances (Table IV). The model efficiency measures show that the model performed best in 2008 and had least skill in 2007. However, the relationship between the observed and mean simulated DOC for all 5 years shows a strong correlation (r 2 = 0.65, Figure 10). Simulations tend to underestimate DOC at higher concentrations, corresponding with DOC exports during the spring melt periods. The results of the sensitivity analysis (Table V) showed that K loss was the most sensitivity model coefficient for the simulations while the distribution of the coefficient was quite narrow. In contrast, the production coefficient (K prod ) was not as important (ranked 6).

DISCUSSION AND WIDER IMPLICATIONS
Understanding the main controlling processes on stream discharge and DOC exports from sub-arctic catchments using a coupled hydrological-biogeochemical model High-latitude regions are highly vulnerable to climate change and are predicted to lose significant proportions of permafrost. It is widely understood that permafrost, seasonal soil frost and slope aspect have a large influence on the hydrology in these catchments (Quinton and Marsh, 1999;Carey and Woo, 2001;Quinton et al., 2005). In addition, biogeochemical processes controlling the production and export of DOC are closely linked to the hydrological processes of the catchment (Carey, 2003;Lyon et al., 2010;Koch et al., 2013). High-latitude catchments experience strong seasonality, and this seasonality is present in the processes and dynamics of the production and export of DOC, where soil moisture and temperature are the major controls.
In this study, a coupled hydrological-biogeochemical model has been developed to simulate stream discharge and DOC exports in a sub-arctic headwater catchment, based on the perceived dominant controlling processes while minimizing model complexity. The model ( Figure 2) presents a simple representation of the dominant hydrological and DOC processes in the catchment. The model builds on the standard HBV model (Lindström et al., 1997) by including components to represent the effects of permafrost, seasonal soil frost and aspect on the hydrological dynamics. In addition, the model incorporates DOC production and loss processes, which are controlled by functions of soil moisture and air temperature.
By including slope aspect and soil freeze-thaw processes, the hydrological component of the model was able to adequately simulate spring melt event duration and magnitude, inter-annual differences of event magnitude and dynamics in late summer, and early autumn and low flows during the winter months. Furthermore, the dynamics of the simulated DOC matched observed DOC concentrations quite well for most study years (r 2 = 0.65 for all years). Despite difficulties in simulating the timing and the magnitude of some of the spring melt events, the observations of each year were within the simulated uncertainty bounds (Figure 3).
During winter, stream discharge is very low because of the presence of permafrost, seasonal frost and a reduction in overall storage. As snowmelt begins, melt water infiltrates into highly permeable organic soils and transports the soil water DOC to the stream. Runoff during the snowmelt event is generally the largest observed throughout the year and leads to the highest stream water DOC concentrations (Carey, 2003;Koch et al., 2013). As the summer season progresses, soil thaw increases and allows for the infiltration of water into the deeper mineral soils, which reduces the soil moisture in the upper soils associated with DOC production (Lyon et al., 2010;Koch et al., 2013). However, model simulations showed that the system remains highly responsive to rainfall events as storage is limited because of the presence of permafrost, which leads to increases in DOC stream concentrations.
In many high-latitude and alpine energy limited systems, the timing of hydrological and biogeochemical processes are dependent on the hillslope aspect, with south-facing slopes experiencing snowmelt and soil thaw before north-facing slopes (Boyer et al., 2000;Dornes et al., 2008;Kicklighter et al., 2013). The model developed here was useful to investigate the role of aspect on DOC production. The simulated DOC in soil water at north-facing and south-facing slopes closely resembled the observations in 2003 (Figure 8), with a rapid reduction in DOC associated with the spring melt, supporting previous research indicating that soil moisture and air temperature as dominant controls of DOC production (Michalzik et al., 2003).

Development and calibration of a simple, coupled hydrological and biogeochemical model
Studies modelling stream water DOC commonly calibrate model parameters of hydrological and biogeochemical model components separately to identify the best hydrological and biogeochemical model parameter sets (Michalzik et al., 2003;Futter et al., 2007;Jutras et al., 2011). Here, we used a multi-objective function including hydrological and biogeochemical model efficiencies simultaneously. The aim of this approach was not to identify a single optimal solution but find a range of plausible parameters for each of the two model components Dick et al., 2014). Stream discharge in regions influenced by snowmelt and permafrost varies greatly throughout the year, with very low flows during winter and very large flows occurring over short periods during the spring melt event, making calibration a challenge (Hamilton, 2001). Using a combination of model efficiency other than the commonly used Nash-Sutcliffe measure, which is biassed towards peak observations, Criss and Winston (2008) and Gupta et al. (2009) allowed for identification of a suitable parameter range for the hydrological component, which does not bias the results of model calibration towards high flow periods.
One of the largest limitations of model development and calibration for hydrological and biogeochemical processes in high-latitude headwater catchments is the lack of empirical observations . Data in these regions is limited because of the difficulties in the 5392 J. S. LESSELS ET AL. collection of observations in these remote areas (Bishop et al., 2008;Holmes et al., 2008) including technical, logistical and financial constraints. Many studies have focussed on the spring melt period because the majority of stream discharge and DOC is exported during this short time (Carey, 2003;McClelland et al., 2007). In our study, DOC observations were restricted to early spring periods associated with the freshet. We acknowledge that by only having observations of DOC during the spring melt, the calibrated parameter range of the model may not reflect processes operating during late summer and autumn, limiting its applicability during this time. However, the ability of the model to reflect the processes during spring and early summer does support the model structure.
The behavioural parameter sets of the coupled model (Table IV) showed a large range of model efficiency measures for DOC and stream discharge individually, which are associated with the large uncertainty bounds during the spring melt period. In addition, the sensitivity analysis revealed that the K loss coefficient was the most sensitive. The K loss coefficient is used to restrict the accumulation of DOC in storage (Equation (6)). As the catchment undergoes substantial accumulation of DOC over the winter period, it is important that the conversion of DOC to other forms of carbon is correctly balanced with the production of DOC, and this assumption is shown by the sensitivity of this coefficient. However, the model was able to capture the main dynamics of both DOC and stream discharge throughout the study period.
Simulated stream discharge reflected the strong snowmelt influence observed during each spring and simulated stream DOC peaked during the stream discharge rising limb (e.g. Figure 9), which has been observed in this system (Carey, 2003). This strong influence of the snowmelt on the system is also reflected in the sensitivity analysis with the SFCF being the second most-important coefficient. The sensitivity of the SFCF is expected as this parameter was used to correct for the potential underestimation of observed snowfall. These results highlight the difficulty of calibrating hydrological and biogeochemical processes in these systems where extreme events (spring melt) have a large influence on the measures of goodness of fit used in the calibration process.

Evaluation of coupled models for representing dominant hydrological and biogeochemical processes
Modelling provides a framework to help overcome the lack of empirical information and to test hypotheses regarding the dominant controls of hydrology and biogeochemistry in these systems. With limited data to inform the modelling, it is important that the output and complexity of the model is appropriate to the data available (Beven and Freer, 2001;Beven, 2006). The coupled parsimonious model developed here captures the main dynamics seen in the system quite well with a r 2 = 0.65 for all DOC observations and a maximum VE of 0.65 for a simulated realization of discharge within each year. Importantly, the model captures prediction of DOC dynamics prior to hydrological events, which is beyond the skill of simple regression based models.
Our model provides a potential framework that could be used to examine the possible effects of climate change in these regions. There are, however, limitations as results suggest that the model struggles to capture the precise timing and magnitude of the spring melt event. In addition, the simulations of both DOC and stream discharge are often associated with large uncertainties. More research is required to investigate the controls of these processes at different scales, focussing on the processes controlling both DOC and stream discharge in these regions. Furthermore, DOC observations are required during late summer to confirm the model DOC simulations. These limitations highlight the difficulties of capturing the dominant processes in complex systems; however, this should not deter the use of models to investigate the controls in these systems.

CONCLUSION
The extent of permafrost in high-latitude regions is predicted to decrease due to climate warming (McGuire and Anderson, 2009;Schaefer et al., 2011). The effect of these changes on hydrological and DOC dynamics in headwaters remains unclear, and it is important that the dominant controls are understood to assess the potential impact of these changes (Kawahigashi et al., 2004;McClelland et al., 2007). In remote headwater catchments in such high-latitude environments, the collection of data is often challenging, and most studies must rely on limited data to investigate the dominant biogeochemical and hydrological processes.
We presented a coupled hydrological-biogeochemical model to simulate discharge and DOC dynamics in highlatitude headwaters. Results show that the parsimonious model was able to capture the dominant controls on stream discharge and DOC in this headwater catchment. However, future work is required to help reduce the uncertainty in both hydrological and biogeochemical components. This could be achieved, for example, by developing models that incorporate information gained from using hydrological tracers for this region. In addition, observations should be expanded to capture DOC exports during late season.

5393
A COUPLED HYDROLOGY-BIOGEOCHEMISTRY MODEL