Carbon Budget Estimation of a Subarctic Catchment Using a Dynamic Ecosystem Model at High Spatial Resolution

A large amount of organic carbon is stored in high-latitude soils. A substantial proportion of this carbon stock is vulnerable and may decompose rapidly due to temperature increases that are already greater than the global average. It is therefore crucial to quantify and understand carbon exchange between the atmosphere and subarctic/arctic ecosystems. In this paper, we combine an Arctic-enabled version of the process-based dynamic ecosystem model, LPJ-GUESS (version LPJG-WHyMe-TFM) with comprehensive observations of terrestrial and aquatic carbon fluxes to simulate long-term carbon exchange in a subarctic catchment at 50 m resolution. Integrating the observed carbon fluxes from aquatic systems with the modeled terrestrial carbon fluxes across the whole catchment, we estimate that the area is a carbon sink at present and will become an even stronger carbon sink by 2080, which is mainly a result of a projected densification of birch forest and its encroachment into tun-dra heath. However, the magnitudes of the modeled sinks are very dependent on future atmospheric CO 2 concentrations. Furthermore, comparisons of global warming potentials between two simulations with and without CO 2 increase since 1960 reveal that the increased methane emission from the peatland could double the warming effects of the whole catchment by 2080 in the absence of CO 2 fertilization of the vegetation. This is the first process-based model study of the temporal evolution of a catchment-level carbon budget at high spatial resolution, including both terrestrial and aquatic carbon. Though this study also highlights some limitations in modeling subarctic ecosystem responses to climate change, such as aquatic system flux dynamics, nutrient limitation , herbivory and other disturbances, and peatland expansion , our study provides one process-based approach to resolve the complexity of carbon cycling in subarctic ecosystems while simultaneously pointing out the key model developments for capturing complex subarctic processes.


Introduction
The high latitudes are experiencing greater temperature increases than the global 5 average (AMAP, 2012;IPCC, 2013). Low decomposition rates due to the cold environment have led to an accumulation of large carbon (C) pools in litter, soils and peatlands, much of which is currently held in permafrost (Tarnocai et al., 2009;Koven et al., 2011;Hartley et al., 2012). However, these C stores may be mineralized rapidly to the atmosphere due to the warming effects on soil microbial 10 activity and thereby increase atmospheric concentrations of both carbon dioxide (CO 2 ) and methane (CH 4 ). Meanwhile, temperature-induced vegetation changes may mitigate those effects by photosynthetic enhancement, which is, however, greatly influenced by disturbances such as plant-herbivore interactions as well as soil water and nutrient contents (Jonasson and Michelsen, 1996;Van Bogaert et al., 2009;15 Keuper et al., 2012). It is becoming crucial to understand those aspects of vulnerable high latitude ecosystems and their responses to climate warming (Callaghan et al., 2013). Ecosystems fix atmospheric C through photosynthesis and return this C back through diverse paths and in different forms. In recent years, many field measurements have been carried out in subarctic and arctic environments to quantify C exchanges 20 between the atmosphere and the biosphere. Those measurements enable us to better understand possible feedbacks from terrestrial biota and responses to the changing climate (Bäckstrand et al., 2010;. However, some concerns about field measurements in the subarctic/arctic environment have been raised and the following research needs have been identified: 1. Complete year-around observations are generally missing and many studies focus only on growing season measurements (Grogan and Jonasson, 2006;936 Roulet et al., 2007;McGuire et al., 2012). Year-around observations are needed because there is clear evidence that C fluxes in the cold seasons are very important (Larsen et al., 2007b;Mastepanov et al., 2008).
2. Observations of interactions between terrestrial and aquatic systems are lacking (Lundin et al., 2013;Olefeldt et al., 2013). Quantifications of terrestrial lateral loss 5 of C are needed not only because they represent a significant fraction of net ecosystem exchange but also because they are intrinsically linked to downstream aquatic C cycling (Lundin et al., 2013). Integrating terrestrial and aquatic C cycling is of high importance for our understanding of the C balance at the catchment scale, particularly at high latitudes. Northern peatlands are large sources of Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and highly instrumented area, can be applied to vast areas where large C stocks exist but long term measurements are lacking. The Stordalen catchment contains several distinct landcover types, including tundra heath, birch forest with heath understory, peatlands and lakes/rivers. Hereafter, the peatland is divided into three groups named "palsa", "Sphagnum site" and "Eriophorum site", based on surface hydrology, permafrost condition as well as characteristic plant communities. An earlier compilation of the C balance of the larger Torneträsk catchment, which encompasses the Stordalen catchment, indicated that there was a significant sink capacity in the birch forest as well as across the peatland (Christensen et al., 2007). This assessment, however, lacked year-around measurement of CO 2 10 and CH 4 emissions and did not include direct measurements of aquatic C fluxes (Christensen et al., 2007). Subsequent observations in the Stordalen catchment have focused on filling the missing components. Consequently, recently updated year-around CO 2 and CH 4 measurements in the peatland identified the wetter nonpermafrost Eriophorum sites to be strong C sinks (−46 g C m −2 yr −1 ) with high CH 4 15 emissions (18-22 g C m −2 yr −1 ) (Jackowicz- Korczyński et al., 2010;, while measurements conducted on the drier palsa (where permafrost is present) showed a relatively weaker uptake (−39.44 g C m −2 yr −1 ) . Total waterborne C exports (DOC plus particulate organic C (POC) and dissolved inorganic C (DIC)) from the terrestrial ecosystems (both peatland and forest) were also monitored 20  and found to represent a significant component of the net ecosystem C balance, ranging from 2.77 to 7.31 g C m −2 yr −1 . In contrast, four years of continuous Eddy-Covariance (EC)-tower based CO 2 measurements in the birch forest revealed very variable C sink/source functionality, which in two out of four years has been found to be a C source to the atmosphere (Heliasz, 2012). Tundra heath, 25 another important landcover type in this region, has lower C uptake (Christensen et al., 2007;Fox et al., 2008) and both birch forest and tundra heath were found to have high spatial heterogeneity (Fox et al., 2008;Heliasz, 2012). Altogether, different landcover types show diverse contributions to this subarctic ecosystem's C balance. With pronounced future warming expected in this region, the structure and function of the different vegetation types are expected to vary dramatically as has been observed during warming in the recent past (Callaghan et al., 2013). In addition, changes to soil conditions due to warming and permafrost thaw will likely stimulate C fluxes to the atmosphere and affect the long-term accumulated C (Wolf et al., 2008a;McGuire et al., 5 2012), but this likely C release needs to be weighed against the possibility of increased uptake by increased primary productivity resulting from longer growing seasons and/or potential CO 2 fertilization. The field measurements described above provide an insight into the ongoing processes and current ecosystem status, but until now, no modeling exercises have 10 been implemented in this region in combination with the comprehensive measured data. Moreover, high spatial resolution predictions of future potential dynamics of both vegetation and soil processes and their responses to the projected climate are lacking in this region. In this study, therefore, we aim to assess the Stordalen catchment C budget in a retrospective as well as in a prognostic way by implementing a process-15 based dynamic ecosystem model (Smith et al., 2001;Miller and Smith, 2012) integrated with a distributed hydrology model (Tang et al., 2014a, b) at high spatial resolution (50 m by 50 m, (Yang et al., 2012)). We quantify the overall C budget of the study catchment by synthesizing diverse C fluxes and specifically address the following questions: (1) will this subarctic catchment become a C source or a larger C sink in the 20 near future? (2) How differently will the catchment's vegetation micro-types respond to the climate drivers? And (3) what are the major limitations in the model's prognostic ability?
To answer these questions, we implemented an arctic-enabled version of the dynamic ecosystem model, LPJ-GUESS WHyMe (Wania et al., 2009a(Wania et al., , b, 2010). This 25 model has been widely and successfully implemented for estimating and predicting ecosystem function in high-latitude regions (McGuire et al., 2012;Miller and Smith, 2012;Zhang et al., 2013). LPJ-GUESS WHyMe includes comprehensive process descriptions to capture the interactions between atmosphere-vegetation-soil domains  (Pilesjö and Hasan, 2014). The study presented here is the first modeling exercise to combine all the available year-around measured data in the Stordalen catchment and at high spatial resolution.
2 Model descriptions 10 The process-based dynamic ecosystem model, LPJG-WHyMe-TFM, was chosen as the platform for studying the subarctic catchment C balance. The processes in the model include vegetation growth, establishment and mortality, disturbance, competition between plant individuals for light and soil water (Smith et al., 2001(Smith et al., , 2014 and soil biogeochemical processes (Sitch et al., 2003). These processes are operated in 15 a number of independent and replicate patches. Vegetation in the model is defined and grouped by plant function types (PFTs), which are based on plant phenological and physiognomical features combined with bioclimatic limits (Hickler et al., 2004;Wramneby et al., 2008). Bioclimatic conditions determine which PFTs can potentially grow in study regions, and vertical stand structure together with soil water availability 20 further influence PFTs establishment based on PFTs' shade tolerance and drought tolerance characteristics. A list of the simulated PFTs in this catchment can be found in  Wania et al. (2009b). The model is driven by climate data and includes both non-peatland and peatland hydrological processes (Fig. 1). The vertical water movement between atmosphere, vegetation and soil is based on Gerten et al. (2004) while the lateral water movement 5 between grid cells was implemented by Tang et al. (2014b) based on topographical variations. More recently, an advanced multiple flow algorithm TFM (Pilesjö and Hasan, 2014) was chosen to distribute water among grid cells, due to its better treatment of flow continuity and flow estimation over flat surfaces (Tang et al., 2014a). Within the catchment boundary, surface and subsurface runoff can move from one upslope 10 cell to multiple downslope cells, which greatly improves hydrological flux estimations (Tang et al., 2014a) and results in a better estimate of C fluxes in peatland region. The soil temperature estimation is driven by surface air temperature, and the Crank-Nicolson heat diffusion algorithm (Crank and Nicolson, 1996) is implemented to calculate the soil temperature profile on a daily time-step (Wania et al., 2009a). The

15
C cycling descriptions in LPJG-WHyMe-TFM for peatland cells are based on Wania's developments (Wania et al., 2009a(Wania et al., , b, 2010, while the non-peatland grid cell C cycling is kept the same as in LPJ-GUESS (Smith et al., 2001;Sitch et al., 2003). The full hydrological processes and peatland C cycling descriptions in LPJG-WHyMe-TFM can be found in Sect. S1 in the Supplement. A brief summary of relevant C cycling 20 processes in the model will be presented below.
A modified Farquhar photosynthesis scheme (Haxeltine and Prentice, 1996;Haxeltine et al., 1996;Sitch et al., 2003) is used to estimate gross primary production (GPP), which is related to air temperature (T ), atmospheric CO 2 concentration, absorbed photosynthetically active radiation (PAR) and stomatal conductance. Part of 25 the GPP is respired to the atmosphere by maintenance and growth respiration (R a ), and the remaining part is net primary production (NPP) for each PFT. The reproduction costs are subtracted from NPP, and thereafter the remaining NPP is allocated to different living tissues in accordance with a set of PFT-specific allometric relationships  (Smith et al., 2001). Leaves, fine root biomass and root exudates are transferred to the litter pool with a given turnover rate, and above-ground plant materials can also provide inputs to the litter pool due to stochastic natural disturbance events and mortality (Smith et al., 2001;Thonicke et al., 2001). The majority (70 %) of the litter is respired as CO 2 directly to atmosphere, with a fixed fraction entering into fast-and slow-turnover 5 soil organic pools (fSOM and sSOM) (Sitch et al., 2003). The overall decomposition rate in the model is strongly influenced by soil temperature and moisture (Sitch et al., 2003). Additionally, the model also estimates emissions of biogenic volatile organic compounds (BVOC) per PFT (Arneth et al., 2007;Schurgers et al., 2009). However, the modelled BVOC values are not compared in the current study, because they 10 only represent a very small fraction of the modelled NEE and, additionally, there is insufficient growing season measurement data in the study domain to evaluate model performance. The major C cycling pathways can be found in the solid line-box area of Fig. 1. For peatland cells, an extra potential C pool for methanogens (C CH 4 ) has been added (see dashed line-box area in Fig. 1) and mainly includes root exudates and 15 easily-decomposed materials (Wania et al., 2010). The majority of C CH 4 is located in the acrotelm layer (Sect. S1) and the oxidation and production of CH 4 together determine the net emission of CH 4 . In the model, the oxidation of CH 4 through methanotrophic bacteria is turned into CO 2 , whereas the un-oxidized CH 4 can be released to the atmosphere by plant-transport, diffusion and ebullition (see the lines with a solid arrow 20 in Fig. 1). The oxidation level is mainly determined by the location of water table position (WTP) in the model (Wania et al., 2010). The biomass production in the current version of LPJ-GUESS WHyMe has no representation of nitrogen (N) limitation and neither N fluxes nor C-N interactions are included (Sitch et al., 2007). The latest version of LPJ-GUESS does include N cycling 25 and N limitation on plant production (Smith et al., 2014), but this capability has not yet been integrated with the customized arctic version of the model adopted in the present paper. Moreover, processes determining concentrations of DOC and DIC in soil water have not yet been explicitly described in the model (Fig. 1) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dissolved C losses and gains in our assessment of the catchment C budget, DOC is estimated by combining modelled runoff with observed DOC concentrations. Measured DIC export data is used directly, based on observations by Olefeldt et al. (2013).

Study site and materials
3.1 Stordalen catchment 5 The Stordalen catchment, located in the discontinuous permafrost region, is situated around 9.5 km east of the Abisko Naturvetenskapliga Station (ANS, Abisko Research Station). The catchment, covering around 16 km 2 , has a steep topography on the southern part and turns into the lower flat peatland region to the north (Fig. 2b). The mean air temperature obtained from ANS records was −0.7 • C for the period 1913- in the forest (see Callaghan et al. (2013) for details of vegetation distribution and changes). Around 5 % of the catchment is covered by lakes or small rivers (Lundin et al., 2013). A detailed description of catchment hydrological conditions can be found in Olefeldt et al. (2013). The map in Fig. 2a is based on an object-based vegetation analysis in this area (Lundin et al., 2014). The classification is based on an aerial 5 imagery dataset collected in August 2008, with a spatial resolution of 0.08 m by 0.08 m.

Model inputs
Monthly climate data at 50 m resolution, including temperature, precipitation, and cloudiness, together with annual CO 2 concentration were used to drive the model from climate at the beginning of the study period by using the first 30 years of climate forcing data to repeatedly drive the model for 300 years. To run the model into the future, the monthly anomalies of simulation outputs from the Rossby Centre Atmosphere Ocean (RCAO) regional climate model  for the grid cell nearest Stordalen were estimated and applied to the historical datasets at 50 m, assuming the same 25 anomaly for all cells in the study region. RCAO climate data was downscaled for an arctic domain using boundary forcing from a general circulation model forced by emissions from the SRES A1B scenario for the period 2013-2080 (Zhang et al., 2013). A soil map provided by the Geological Survey of Sweden was used to identify the peatland fraction of each cell. Notably, the detailed rock area shown in Fig. 2a is not represented by the soil map used. Imagery Digital Elevation Model (DEM) from the 5 National Land Survey of Sweden was used and the TFM algorithm (written in a Matlab script, R2012) was used to calculate topographic indices (Pilesjö and Hasan, 2014) used in the model LPJG-WHyMe-TFM.

Additional model parameterization and sensitivity testing
An artificially-adjusted snow density (100 kg m −3 , following the measurements ranges provided in Judson and Doesken, 2000) was implemented for birch forests to mimic deeper trapped snow in the forest. The observed thick top organic soil in the birch forest (Heliasz, 2012) is also represented in the model by increasing the organic fraction of the top 0.1 m soil layer (Olefeldt and Roulet, 2014) in the model. Additionally, the model simulation uses a single replicate patch in each 50 m grid cell. Finally, an additional 15 simulation was performed in which we kept the CO 2 concentration constant after 1960 and allow the remaining inputs to vary as normal to diagnose the extent to which the CO 2 concentration influences on ecosystem dynamics and C fluxes.

Model evaluation data
The observations used to evaluate the model include: (1)    , (6) EC-tower measured birch forest CO 2 NEE for four continuous years, 2007-2010(Heliasz, 2012), (7) catchment-level annual CO 2 and CH 4 emission from lakes and streams from 2008 to 2011 (Lundin et al., 2013). 5 The catchment level C budget was calculated using estimated C fluxes (CF i ) (after first accounting for observed DOC or DIC fluxes, where available) weighted by landcover fractions (f i in Eq. 1). The symbol i in Eq.

Calculation of the catchment carbon budget
(1) represents different landcover types, in our case, birch forest, peatland (Sphagnum and Eriophorum sites), tundra heath and lake/river. To identify the temporal variations of C fluxes, the identification different 10 landcover (except peatland region) is based on the dominant PFTs for each grid cell during the reference period of 2001-2012. The aerial photo-based classification of different peatland types is used when estimating peatland fluxes as a whole .

Results
The modelled C fluxes for birch forest, peatland and tundra heath are first compared with the measured data for the period of observation. The seasonality and magnitudes of the fluxes are evaluated and then used to estimate the catchment-level C budget. Both the plant communities and hydrological conditions in the peatland differ among the palsa, Sphagnum and Eriophorum sites. Notably, the measured NEE (see grey bars in Fig. 3a and b) covers both dry hummock (palsa) and semi-wet (Sphagnum 5 site) vegetation , whereas our model cannot represent the dry conditions of the palsa. Therefore, the observed fluxes over both palsa and Sphagnum site were compared with the modelled fluxes at the Sphagnum site. The C fluxes magnitudes (including both wintertime emission and summertime uptake) are larger at the Eriophorum site, when compared with the palsa/Sphagnum site for both measured 10 and modelled data. The modelled NEE at both sites generally captures the seasonality and magnitude of measured NEE, from being a strong sink (negative NEE) of CO 2 during the summer (mainly June-August) to being a wintertime CO 2 source. However, the model is unable to fully capture the C source/sink functionality in September at both sites. Furthermore, the modelled winter respiration at the Eriophorum site . For the Eriophorum site, the three-year mean growing season uptake of C is underestimated by the model (Fig. 3d), which indicates that the modelled photosynthetic rates may be too low, or summer respiration rates may be too high, or both. 25 Two full years (2006 and 2007) of EC-tower measured CH 4 emissions were used to evaluate modelled estimations, and three different pathways of modelled CH 4 are also presented ( Fig. 3e and f) (20.23 and 22.57 g C m −2 yr −1 for the measured and modelled, respectively). Specifically, the wintertime emissions are slightly underestimated but this underestimation is compensated by the overestimated summer emissions. Plant-mediated transport of CH 4 dominates during the growing 5 season, while the ebullition and diffusion transport reach a maximum in August. Additionally, the plant-mediated CH 4 emission is the main pathway active during the late spring and early autumn.

Birch forest
The modelled average LAI for the birch forest and understory vegetation is around 10 1.4 and 0.3, respectively, which are consistent with observations made in this area (Heliasz, 2012). Modelled and measured monthly and cumulative NEE are compared for the years 2007-2010 in Fig. 4. From the monthly NEE comparisons (Fig. 4a), we see that the model underestimates ecosystem respired C both before and after the growing season. The maximum photosynthesis-fixed C in July is lower than that measured for

Tundra heath
Around 29.8 % of the catchment with alpine terrain was covered by heaths and dwarf shrubs during the reference period 2001-2012. The model predicts that low-growing 10 evergreen shrubs (PFT: LSE, e.g. Vaccinium vitis-idaea) currently dominate in this area, with tall summer-green shrub (PFT: HSS, e.g. Salix spp.) dominant in the future predictions (2051-2080) (Fig. 5). Since there is no available year-around observations of the C balance in this part of the Stordalen catchment, a synthesis of published data from similar environments is used to evaluate the model estimations (Table 1). 15 Four periods with averaged C uptake values from our model are presented (see the last four columns in Table 1) and a clear increase of summer uptake can be found with increased temperature and CO 2 concentration. The modelled, whole-year uptake also follows the increasing trend, except for the period of 2051-2080. The model-estimated summer C uptake is much stronger than the observations made at the high arctic 20 heath of NE Greenland, which is reasonable when considering the longer growing season in our catchment. The study conducted in an alpine area in Abisko by Fox et al. (2008) shows a wide range of estimated NEE over 40 days during the summer time of 2004. The wide range of NEE values are because three levels of vegetation coverage were studied and the lowest uptake is from sparse vegetation dominated 25 by bare rock and gravels, which is similar to the situation in our tundra heath sites (Fig. 2a) 2.33 g C m −2 yr −1 from the observations and 3.09, 2.03, and 1.69 g C m −2 yr −1 from the model estimations, respectively. The downward trend over the three years is captured and is related to lower precipitation during the latter years. The observed DOC export rates at the palsa/Sphagnum and Eriophorum sites in 2008 are directly used to represent different types of peatland export level and the values are 3.35 15 and 7.55 g C m −2 yr −1 for the palsa/Sphagnum and Eriophorum sites, respectively . DIC export is currently beyond the scope of our model but nonetheless contributes to the whole C budget. We used an averaged DIC export value of 1.22 g C m −2 yr −1 based on the published data in Olefeldt et al. (2013) and DIC export is included in the estimation of the catchment C budget below. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | CH 4 and CO 2 , but the streams dominate CO 2 emission while the lakes dominate CH 4 emissions. Averaging across the catchment area, the measured annual CO 2 and CH 4 emissions from streams are 10.1 and 0.06 g C m −2 yr −1 , respectively, while the lake emitted 0.5 g C m −2 yr −1 as CO 2 and 0.1 g C m −2 yr −1 as CH 4 (Lundin et al., submitted).
Since river and lake processes are at present beyond the scope of the model, the 5 contribution of aquatic systems to the catchment C fluxes is purely based on the observed data.

Modelled whole-catchment carbon balance
One of the benefits of using a dynamic vegetation model is that it allows us to investigate the vegetation and C budget responses to climate change. In this section, 10 annual variations of both climate variables (temperature and precipitation) and different ecosystem C fluxes are presented for the whole simulated period. In addition, the normalized catchment C budget is also estimated.
Our simulation results suggest that the temperature increase (around 2 • C, see  Fig. 6-6a-c). The DOC export from birch sites also increases slightly over the study period ( Fig. 6-3a-c). Besides leading to an increase in ecosystem biomass in grid cells currently occupied by birch forest (Fig. 5c and d), warming will also result in the current tundra heath close to the birch treeline being replaced by an 20 upward expansion of the birch forest ( Fig. 5a and b). This is indicated by the changes to birch treeline's uppermost elevation in Fig. 6-8a-c. The dramatic increase of C uptake in the tundra heath region (Fig. 6-5a-c) is largely a result of this vegetation succession. Meanwhile, the successive degradation of permafrost and slightly higher annual precipitation may result in more anaerobic conditions for the modelled peatland. 25 Combined with warmer soil conditions, these result in both higher decomposition rates of soil organic matter and greatly increased CH 4 production. Moreover, a stronger response of respiration than NPP to the temperature increase reduces the net C 951 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | uptake in the peatland (see Fig. 6-4a-c). The catchment, as a whole, shows an increased C uptake (see the trend line in Fig. 6-7a) over the 1913-2080 period. The averaged C budgets for the two selected periods (1961-1990 and 2051-2080) are −1.47 and −11.23 g C m −2 yr −1 , respectively, with the increase being dominated by the large increase seen in the uptake rates at the birch forest and tundra heath sites. 5 During the 2051-2080 period, both boreal needle-leaf shade-tolerant spruce (PFT: BNE, e.g. Picea abies) and boreal needle-leaf (but less shade-tolerant) pine (PFT: BINE, e.g. Pinus sylvestris) start to appear in the birch-dominated regions (Fig. 5). In the current tundra regions, the coverage of HSS greatly increases at the expense of LSE. For the northern parts of the catchment, birch forest densification is observed with future warming, while the greatest relative changes in biomass occur near the treeline. The increased temperature together with increased CO 2 concentration by 2080 are very likely to increase CO 2 uptake in both tundra heath and forested areas, though the nutrient limitations are not included in this version of the model (but see Smith et al., 2014). A summary of the modelled catchment C components and C fluxes from 15 different sources is given in Table 2. Additionally, the modelled estimations during the warm period of 2000-2005 are also presented in Table 2 and the positive value of catchment C-budget is seen to be an exception compared to the main reference period.
To illustrate the impact of CO 2 increases alone on the modelled C budget, the differences between simulations (∆C fluxes) with and without a CO 2 increase since 20 1960 are shown in Fig. 7. The original outputs of those two simulations are plotted separated in Fig. S1 in the Supplement. Interestingly, the simulation with constant CO 2 forcing after 1960 greatly reduces the birch and tundra uptake (positive values in Fig. 7), whereas the peatland NEE and CH 4 are not strongly influenced by the CO 2 increase. The catchment C budget dynamics are consistent with the changes seen in the birch 25 and tundra heath regions.
To account for the fact that CH 4 is a more potent greenhouse gas than CO 2 , an estimation of the global warming potential (GWP) of the two simulations can be made. Assuming the relative climate impact of CH 4 is 28 times greater than CO 2 over a 100 year period (IPCC, 2013), the calculated GWP for the simulation with atmospheric CO 2 increase is three times larger in the period 1961-1990 (27.3 g CO 2 -eq) than the period 2051-2080 (8.8 g CO 2 -eq). However, in the simulation without CO 2 increase, the GWP in 1961-1990 (40.7 g CO 2 -eq) is approximately half of the GWP value for the 2051-2080 period (93.3 g CO 2 -eq). This shows that the change in global warming potential 5 found in the model simulations is strongly influenced by the CO 2 concentration used to force the model, as the CO 2 trajectory can alter the balance between the GWP changes resulting from C uptake in the birch forest and tundra, and the peatland CH 4 emissions.

10
To our knowledge, our model simulation is the first attempt to create a C budget of a subarctic catchment based on a dynamic global vegetation model applied at the local scale and to predict its evolution in response to changing climate. The C budget estimations in this paper include major flux components (CO 2 , CH 4 and hydrological C fluxes) based on a process-based approach. The magnitude and seasonality of 15 modelled C fluxes compare well with the measurements (Figs. 3 and 4 and Tables 1 and 2), which gives us confidence in the ability of our model to represent the main processes influencing the C balance in this region. Our hope is that, by using a processbased modelling approach at high spatial resolution, our methodology will be more robust in estimating the C budget in a changing future climate than other budget 20 estimation methods (Christensen et al., 2007;Worrall et al., 2007). In response to a climate warming scenario, our model shows a general increase in the C sink in both birch forest and tundra heath ecosystems, along with greater CH 4 emissions in the catchment (Fig. 6). Integrated over the catchment, our modelled C budget indicates that the region will be a greater sink of C by 2080, though these estimates are sensitive to the atmospheric CO 2 concentration. Nevertheless, the current model set-up and simulations still contain limitations in both the historical estimations and predictions. Below we discuss the model's current performance as well as a few existing limitations, and further propose some potential model developments. Most of the studies referred to below have been conducted specifically in the area near the Stordalen catchment or in other regions with similar environments.

Peatland CO 2 and CH 4 emission
The model estimations in the peatland demonstrated skill in representing the relative differences of CO 2 fluxes between the palsa/Sphagnum site and Eriophorum site. The current model version cannot represent the very dry conditions of the palsa due to the acrotelm-catotelm soil structure (Wania et al., 2009a); therefore, the 10 overestimated emissions rates of CO 2 ( Fig. 3a and b) are expected when using the modelled fluxes from the Sphagnum site to compare with the measured data covering both the palsa and Sphagnum sites (Larsen et al., 2007a;Bäckstrand et al., 2010). The cold-season respiration at the Sphagnum site can account for at least 22 % of the annual CO 2 emissions (Larsen et al., 2007a), which is a further reason 15 for the model overestimations. Moreover, at the Eriophorum site, the model slightly underestimated summertime uptake, which was found to be related to modelled summer respiration being higher than dark chamber-measured respiration, though with high spatial variations (Tang et al., 2014a). The accurate representations of annual CH 4 emissions at the Eriophorum sites 20 reflect the model's improved estimations of both hydrological conditions and dynamic C inputs to the C CH 4 pools (Fig. 1). With the inclusion of the distributed hydrology at 50 m resolution (Tang et al., 2014a, b), the lower-lying peatland can receive water from the upslope regions and the fact that water can rise to a depth of 10 cm above the surface both creates anaerobic conditions in the model and favours the establishment  (Wania et al., 2010). Nevertheless, it remains the case that excluding the palsa type could result in a general overestimation of CH 4 emission over the whole peatland. 5 The high wintertime CO 2 fluxes to the atmosphere observed in the birch forest and the underestimations by our model (Fig. 4) highlight the importance of the representation of winter season C fluxes for the birch forest, particularly as winter temperatures (T ) are increasing more than summer T (Callaghan et al., 2010). In the model, the respiration rate is linked to soil T and moisture and the underestimated respiration could be attributed to the lower T estimations in the wintertime (covering all months except June-September) when compared with the observed soil temperature at 5 and 10 cm depth in 2009 and 2010 (the modelled wintertime soil T is around 1.9 • C lower than that observed at both depths). Also, the investigations of snow depth insulation influence on soil T and respiratory activity by Grogan and Jonasson (2003) 15 have shown that soil T contributes to the greatest variations in respiratory activity at our birch sites. The birch forest is located on the relatively lower parts of the catchment (Fig. 2) and traps the wind-shifted lighter snow from the upland tundra heath, creating a much thicker snow pack and therefore significantly increasing the soil T (Groendahl et al., 2007;Larsen et al., 2007b;Luus et al., 2013). Snow depth 20 measurements in the Abisko birch and tundra heath sites from 24 March to 7 April 2009 (http://www.nabohome.org/cgi-bin/explore.pl?seq=131) revealed the snowpack in the birch forest was 26.62 cm deeper on average than the snowpack in the tundra heath. The snow density for the birch forest was artificially decreased to 100 kg m −3 in the model in order to increase the snow depth in the absence of a wind-redistribution 25 mechanism in our model, but it is still hard to capture the high soil T as well as high emission rates in winter. The thickness of organic soil layer is another crucial component controlling birch site respiration in our catchment (Sjögersten and Wookey, 2002;Heliasz, 2012). The thickness of organic soil is most likely connected to the past transformation of the sites from heath to forest and is expected to decrease if birch biomass continues to increase (Hartley et al., 2012(Hartley et al., , 2013. Furthermore, Hartley et al. (2012) showed evidence of the 5 decomposition of older soil organic matter during the middle of the growing season and concluded that with more productive forests in the future, the soil stocks of C will become more labile. From a model perspective, faster turnover rates of soil C pools (Sitch et al., 2007) could be used in climate warming model experiments to reflect the accelerated C mobilization in the soil. This could, to some degree, offset the stronger 10 C uptake in the birch forest site in the simulation results shown here. In addition, the current parameterization of the organic fraction in the birch forest may not fully represent the organic horizons in reality soil.

BGD
The birch forest is cyclically influenced by the outbreak of moths (E. autumnata) and an outbreak in 2004 greatly affected the birch forest, resulting in a much lower 15 C uptake than the average (Heliasz et al., 2011). Even though it is assumed that the fluxes during 2007-2010 had returned to the pre-defoliation levels (M. Heliasz, personal communication, 2013), it is still difficult to completely exclude the influence of insect disturbance in this forest. For the current model simulation, stochastic mortality and patch-destroying disturbance events have been included to account for the impacts 20 of these random processes on ecosystem C cycling. However, to give a more accurate and reliable representation of the C fluxes in the birch forest, an explicit representation of insect impacts, outbreaks and their periodicity should be included (Wolf et al., 2008b) as well as other disturbances (Callaghan et al., 2013;Bjerke et al., 2014). For example, a warm event in winter 2007 caused a 26 % reduction in biomass over an area of over

Benefits of high spatial resolution and limitations from monthly temporal resolution
Subarctic ecosystems are characterized by small-scale variations in vegetation composition, hydrological conditions, nutrient characteristics and C fluxes (Lukeno and Billings, 1985;McGuire et al., 2002;Callaghan et al., 2013). The spatial resolution of 5 50 m in this model application allowed us to capture the diverse vegetation micro-types and their altitudinal gradient in the catchment, as well as their differential responses to climate changes (Figs. 5 and 6). This is unlikely to be well represented by a simpleaveraging approach across the landscape, e.g. RCMs or GCMs. Furthermore, it is worth noting that the carbon cycling in the peatland, especially CH 4 fluxes, is sensitive 10 to the WTP estimations (Tang et al., 2014a). Without spatially-distinguished climate and topographical data, it becomes impossible to implement our distributed hydrological scheme and thereby capture the peatland WTP dynamics. To our knowledge, the use of 50 m climate data as forcing of LPJG-WHyMe-TFM (Yang et al., 2011) is among the most detailed and comprehensive modelling exercises related to C cycling. Although 15 this study focuses on C cycling, the innovative projects changes of vegetation microtypes at high spatial resolution are relevant to local stakeholders such as conservation managers and reindeer herders. Nevertheless, the high spatial resolution data is currently produced with coarser monthly temporal resolution, which could restrict the model's ability to accurately 20 estimate C fluxes at the start and end of the growing season (Figs. 3 and 4). Due to the dramatic variations of day length at high latitudes, a few days of misrepresenting the starting date of the growing season could significantly alter the estimated plant C uptake (Heliasz, 2012

Projection uncertainties
A temperature increase of 2 • C (Fig. 6) and elevated CO 2 concentrations could greatly increase vegetation growth and thus the C sink of the whole catchment. However, the densification and expansion of birch forest as well as the increased presence of boreal spruce and pine PFTs in our projected period (Fig. 5) could be 5 strongly influenced by reindeer grazing and herbivore outbreaks (Hedenås et al., 2011;Callaghan et al., 2013), even though those detected changes are consistent with other models simulations (Wolf et al., 2008a;Miller and Smith, 2012) and general historical trends (Barnekow, 1999). Furthermore, climate warming may favor the spread of insect herbivores, so an assessment of ecosystem responses to future climate change cannot 10 ignore these disturbances (Wolf et al., 2008b) and also other factors such as winter warming events (Bjerke et al., 2014) Temperature increase results in a larger extent of permafrost degradation in the future. Meanwhile, the increased amount of available water from precipitation and lateral inflow may increase the degree of anoxia and further favors the flood-tolerant 15 WetGRS PFT growth as well as CH 4 production. However, the exact extent of wetting or drying of peatland in the future is still highly uncertain, and the model prediction depends strongly on the climate scenario, permafrost thawing and the resulting balance between increased water availability and increased evapotranspiration. If the peatland drying is large enough, the reduced degree of anoxia could reduce CH 4 emissions in 20 the future (Wrona et al., 2006;Riley et al., 2011). In our case, the peatland, located at the lower part of the catchment, receives water from the southern mountain region, and is more likely to become wetter in the near future in response to the increases in water availability (Wrona et al., 2006). Additionally, the determination of grid cell peat fractions in our simulation is dependent on the current (historical) soil map; therefore, 25 the projections of peatland expansion to non-peatland cells cannot be reflected in the current model predictions. This could bring some additional uncertainties into the catchment C-budget estimations (Malmer et al., 2005;Marushchak et al., 2013) To cover all the major C components in this catchment, the current estimations of the C budget used available, observed aquatic emissions and DOC concentration and DIC export values, all of which were assumed constant for the whole simulation period, 1913-2080. However, the direction and magnitude of changes to aquatic C fluxes are hard to quantify without modelling additional processes in these systems. 5 Previous studies have found that the substantial thawing of permafrost as well as increased precipitation in recent decades have significantly increased total organic carbon (TOC) concentrations in lakes (Kokfelt et al., 2009;Karlsson et al., 2010). The increased loadings of nutrients and sediments in lakes are very likely to increase productivity of aquatic vegetation, but the effects may be offset by the increased inputs 10 from terrestrial organic C and increases in respiration (Wrona et al., 2006;Karlsson et al., 2010). Emissions of CH 4 from aquatic systems are more likely to increase due to the longer ice-free season (Callaghan et al., 2010), increased methanogenesis in sediments and also increased CH 4 transport by vascular plants (Wrona et al., 2006). Furthermore, a recent study by Wik et al. (2014) found that CH 4 ebullition from lakes 15 is strongly related to heat fluxes into the lakes. Therefore, future changes to energy fluxes together with lateral transports of dissolved C from terrestrial ecosystems to the aquatic ecosystem are especially important for predicting C emissions from aquatic systems.
As we have discussed, the dynamics of birch forest and to a lesser extent tundra 20 heath C assimilation largely determines the catchment's C budget (Figs. 5 and 6), whereas the dramatic increase of CH 4 can slightly offset the net climate impact of the projected C uptake. Furthermore, both modelled C budget and GWP values are very sensitive to the atmospheric CO 2 levels. However, to date, there is no clear evidence showing significant long-term CO 2 fertilization effects in the arctic region 25 (Oechel et al., 1994;Gwynn-Jones et al., 1997;Olsrud et al., 2010). Two possible reasons for this lack of CO 2 fertilization response might be that CO 2 levels close to the moss surface in birch forest can reach as high as 1140 ppm already and are normally within the range of 400-450 ppm (Sonesson et al., 1992) Gehrke, 1996;Johnson et al., 2002), temperature and growing season length (Heath et al., 2005) and forest longevity (Bugmann and Bigler, . Therefore, more field experiments are urgently needed in order to quantify and understand the CO 2 fertilization effects on the various vegetation micro-types of the subarctic environment, and particularly tall vegetation types. Overall, the current model application has been valuable in pointing to this gap in process understanding and meanwhile shows the importance of including vegetation 10 dynamics in studies of C balance. Furthermore, a current inability to include the potential impacts of peatland expansion, potential increases of emissions from aquatic systems as well as the potential nutrient limitations on plants (but see Smith et al., 2014) and disturbances (Bjerke et al., 2014), make it likely that our projections of the catchment C-budget and CO 2 -GWP will vary from those that may be observed in the 15 future. However, our high spatial resolution, process-based modelling in the subarctic catchment provides an insight into the complexity of responses to climate change of a subarctic ecosystem while simultaneously revealing some key uncertainties that ought to be dealt with in future model development. These developments would be aided by certain new observations and environmental manipulations, particularly of Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | for the providing us with regional climate projection data and anomalies. PAM acknowledges financial support from the ADSIMNOR (Advanced Simulation of Arctic Climate and Impact on Northern Regions) project funded by FORMAS (Swedish Research Council), and the Lund University Centre for the study of Climate and Carbon Cycle (LUCCI) funded by VR. The study is a contribution to the strategic research area Modelling the Regional and Global Earth System (MERGE), the Nordic Centre of Excellence DEFROST and the EU PAGE21 project. TVC and TRC wish to thank FORMAS for funding for the project "Climate change, impacts and adaptation in the sub-Arctic: a case study from the northern Swedish mountains" (214-2008-188) and the EU Framework 7 Infrastructure Project "INTERACT" (http://www.eu-interact.org/).  Period 1997-2007, 2007, 13 Jul-1961-1990, 1991, 2051-2080, 2005   Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 6. Time-series (column a) of annual mean temperature, precipitation sum, NEE (g C m −2 yr −1 ), averaged DOC export (g C m −2 yr −1 ) from the birch forest, birch tree line elevation (m) and catchment C-budget (g C m −2 yr −1 ) values between 1913 and 2080. The NEE data for the birch forest and the peatland have had the corresponding DOC and DIC values subtracted.

BGD
Here the averaged C fluxes from the Stordalen peatland (the northeast side of the catchment) are used to represent for the averaged C fluxes for the whole catchment peatlands since only the Stordalen peatland DOC and DIC observations are currently available. The aerial photo-based classification of the Eriophorum and Palsa/Sphagnum peatland fractions within the Stordalen peatland is used to scale up the C fluxes. The trend of each dataset is shown with a red dashed line. The second and third columns (b and c) of the figure focus on the periods 1961-1990 and 2051-2080 (these two periods are also indicated in the first column with shaded area). The numbers in bold type in columns (b) and (c) show the annual average of each quantity for the respective time period. The fractions of peatland, birch forest, tundra and lakes/rivers are 5.7, 57.69, 29.76 and 5 %, respectively. There are around 1.6 % of the catchment area are estimated with the dominance of C3G and HSE which are not included in the above classification. The last row shows the birch tree line maximum elevation changes over the period.