Diurnal , seasonal and long-term behaviours of high Arctic 1 tundra heath ecosystem dynamics inferred from model 2 ensembles constrained by the time-integrated CO 2 fluxes 3 4

9 Ecosystem CO2 fluxes in high Arctic are rather dynamic, as they are sensitive to climatic variability through 10 multiple ecosystem processes, for instance, vegetation and snow dynamics as well as permafrost thawing, operating 11 at different time scales. Uncertainties from both high-frequency measurements and model assumptions challenge 12 model calibration to describe both shortand long-term phenomena related to weather and climate variabilities. In 13 this study, we generated three model ensembles using a Monte-Carlo based uncertainty approach with acceptance 14 criteria for 15 years of eddy covariance CO2 measurements of a high Arctic heath ecosystem based on the time15 integrated CO2 fluxes within the day, the year and the entire period. The temporal distribution of residuals between 16 the model and measurements indicated that the three model ensembles reasonably simulated diurnal, seasonal and 17 long-term behaviours of CO2 fluxes respectively. The inter-annual variation of CO2 fluxes over 15 years showed the 18 current ecosystem is at a transition from being a C sink to a C neutral balance. The long-term behaviour model 19 ensemble simulated a more intensified diurnal C cycle than the short-term behaviour model ensembles. The 20 intensified C cycle was mainly attributed to a faster depletion of the soil C pools. The sensitivities of posterior 21 parameters to the model performance index (coefficient of determination, R) reflected that parameters in the 22 processes of soil water and heat transfer and snow dynamics regulated the short-term behaviour of CO2 fluxes, while 23 parameters in the process of soil decomposition regulated the long-term behaviour of CO2 fluxes. Our results 24 suggest that the development of ecosystem models should diagnose their effectiveness in capturing ecosystem CO2 25 exchange behaviour across different time scales. A clear trade-off may exist when the model is tuned to capture both 26 the shortand long-term variation of CO2 fluxes. To constrain the model with the time-integrated CO2 fluxes is a 27 simple and useful method to reduce the non-explained errors and to identify the crucial link to controlling 28 parameters and processes. 29


Introduction
Northern permafrost regions characterized by cold climate, short growing seasons and poorly-drained soils with low fertility hold more than half of global soil carbon (C) with the top 3 m (Hugelius et al., 2014).Currently, widespread evidence reveals that rapid warming of the Arctic has resulted in permafrost thawing (Grosse et al., 2016).This implies that further warming can potentially release significant amounts of soil C to the atmosphere, exerting a positive feedback to global climate change (MacDougall 2012; Koven et al., 2015;Schuur et al., 2009Schuur et al., , 2015)).On the other hand, enhanced warming and redistributed wetting may alter growing season length, vegetation phenology and nutrient and water availability to conditions favored by shrubby and productive ecosystems, which may sequester more carbon dioxide (CO2) comparing to present tundra ecosystems (Xia et al., 2017;Zhang et al., 2013aZhang et al., , b, 2014;;Elmendorf et al., 2012;Myers-Smith et al., 2015).These two opposing ecosystem feedbacks pose a challenge to identify current and future transition of Arctic ecosystems functioning between C sink and source.
The net C gain or loss in ecosystems is described by net ecosystem exchange (NEE) of CO2 between the atmosphere and biosphere.Its variability is affected by responses of two constituent fluxes (i.e.gross primary production (GPP) and ecosystem respiration (ER)) to abiotic and biotic drivers across a wide range of time scales (hours to years) (Richardson et al., 2007;Baldocchi et al., 2017;Wu et al., 2017).Previous studies used models or wavelet spectra analysis to disentangle the effects of short-and long-term climatic variation on eddy flux variability (Baldocchi et al., 2001;Braswell et al., 2005;Stoy et al., 2005Stoy et al., , 2009Stoy et al., , 2013;;Richardson et al., 2007).They generally agreed that climatic variation tightly controls the short-term variation of CO2 fluxes, and the effects will become progressively less important at a long term.The intermediate scales of weather patterns, such as freezing events, heat waves or rain pulses, may influence plant phenology and microbial activities through their effects on soil temperature and soil moisture (Vargas et al., 2010).Accordingly, the dominant drivers of ecosystem dynamics can shift from physical or physiological controls to biological responses to climatic variation along multiple temporal scales (Richardson et al., 2008;William et al., 2009;Wu et al., 2012a;Fig. 1).This transition possibly becomes remarkable in the Arctic tundra, which has experienced the substantial changes in seasonal climate and thawing depths of permafrost during the past few decades (AMAP, 2011;Myers-Smith et al., 2015;Schuur et al., 2015).
Soil heat and water transfer, snow dynamics, vegetation dynamics, and soil decomposition are major processes to regulate Arctic ecosystems NEE variation across multiple time scales.Soil heat and water fluxes respond to diurnal climatic variation instantaneously.In situ manipulation experiments, laboratory incubations and modelling studies show that increased soil temperatures due to the effects of air warming or more efficient snow insulation are primary abiotic drivers for increased soil effluxes in both growing seasons and cold seasons (Elberling, 2003;Semenchuk et al., 2016;Webb et al., 2016;Li et al., 2014).Soil moisture is associated with water stress of photosynthesis and sometimes may become a limiting factor other than soil temperature by determining the form and magnitude of soil C release to the atmosphere (Blok et al., 2016;Natali et al., 2015).Snow dynamics is another important aspect that may cause changes in soil temperature, soil moisture, growing season length, and quantity and quality of substrate (Lund et al., 2012;Grøndahl et al., 2007).Snow directly impacts on the onset of growing seasons but also has carryover effects on the summer-and autumn-time CO2 exchange through soil moisture limitation (Westergaard-Nielesen et al., 2017).Moreover, increased soil temperature due to a larger winter snow may also affect the nitrogen (N) 2 Materials and methods

Site description
The study site is a Cassiope tetragona heath tundra situated at a lower part of Zackenbergdalen in Northeast Greenland (74°28Ń, 20°34Ẃ, 38 m a.s.l),where the CO2 flux measurements were collected by an EC mast since 2000 (Lund et al., 2012).Within the study area, a climate station was established in summer of 1995 and started to collect the hourly meteorological measurements since then.The local climate type is characterized as "polar desert" with low precipitation in summer and intensive and persistent coldness for the rest part of the year (Abermann et al., 2017).Based on the meteorological measurements, annual mean air temperature shows an increase of 0.05 °C yr -1 for 1996-2014, and such increased warming is largely attributed to the winter-time warming (i.e.DJF) of 0.15 °C yr - 1 .Both changes in annual total precipitation and snowfall are less than 1% during 1996-2014.Annual maximum thawing depth has been observed around 60-70 cm over August and September (Hollesen et al., 2011).C. tetragona, an evergreen dwarf shrub species with the average height of 5-10 cm (Elberling et al., 2008;Campioli et al., 2012) and leaf area index of 0.4-0.5 (Campioli et al., 2009), dominates the surrounding ecosystems by covering around 53% of surface areas (Bay 1998).Mosses and lichens cover the areas around 8% and 3% respectively.Few patches of Salix Arctica, Dryas and other herbs cover less than 5% of the total areas.The surface albedo of summerand winter-time varies from 0.15 to 0.8.The soil type is classified as Typic Psammoturbels, which locates most soil C near the surface and a consistently low C concentration for the deeper soil horizons (Elberling et al., 2008).In some places, the buried organic-rich surface layer at approximately the depth of 20-30 cm was dated to reflect soil development since the Holocene Climate Optimum around 5000 years ago (Christiansen et al., 2002).Table 1 lists the information on climate, plant and soil characteristics, some of which are used to set initial conditions of soil and plant properties in the model and will be mentioned in the following sections.

Measurements and data
The hourly meteorological data sets for 1996-2014 consist of incoming short-wave and long-wave radiation, precipitation, 2m air temperature, wind speed, and relative humidity.The daily soil temperatures were measured at the depth of 0, 60, 130 cm for the same period as the meteorological measurements.Snow depth for 1997-2014 was measured at every three hours.The CO2, sensible heat and latent heat fluxes were measured by the EC technique based on the half-hour average after being gap-filled or storage-corrected for the period 2000-2014 (Lund et al., 2012), but they were integrated into the hourly interval for the model comparison.The 2 depths (10 and 30 cm) of soil water content for 2006-2014 were measured by the Time-Domain Reflectometry (TDR) technique.All these data sets are available in the Greenland Ecosystem Monitoring database (http://data.g-e-m.dk/).

Model description
The CoupModel (version 5.0, available in http://www.coupmodel.com) is a process-oriented ecosystem model with a number of modules that represent biotic and abiotic processes for heat and mass (i.e.water, C and N) transfer within an atmosphere-plant-soil continuum and depth-dependent environmental controls on soil C dynamics (Jansson et al., 2012).The following model description mainly focuses on the processes relevant to this work, and the parameters and equations are provided in the Table S1.A more detailed documentation of principles, parameters and equations can be found in Jansson and Karlberg (2010).

Soil heat and soil water
Soil heat and water flows were estimated using a coupling between the Fourier equations and Richard equations based on soil physical characteristics of a vertical multi-layered domain (Jansson and Halldin, 1979).In this study, the soil profile (37.35 m) was prescribed using 35 soil layers with increasing thickness from 0.05m for the top layer to 3 m for the bottom layer.Soil heat flow only accounted for conduction of heat using the thermal conductivity function.Soil thermal conductivity was estimated by considering soil texture characteristics, organic matter thickness, soil water content, and ice content (Kersten, 1949).The upper boundary conditions for soil heat transfer were regulated by soil surface heat flow, which was calculated using a function of soil surface temperature, organic layer thickness and convective heat flow given by infiltrated water and vapour.The soil surface temperature was a weighted sum of the estimation on three conditions: the bare soil, the soil partially covered by snow and the soil partially covered by plants.This allows soil surface temperature to reflect either the interactions between aerodynamic properties (e.g.plant cover and bare soil) or the steady-state heat transfer between soil and a homogeneous snow pack.The lower boundary condition for heat conduction was estimated by a parameter which represents a constant geothermal heat flow.
Soil water was balanced by water infiltrated into the soil, soil evaporation, plant water uptake and runoff.Water flows between adjacent soil layers were assumed to be laminar, consisting of a matrix flow, a vapour flow and a macro-pores bypass flow.The upper boundary condition for calculating water flows was given by separate subroutines accounting for surface hydrology, soil frost, snow dynamics, and canopies interception of precipitation.The lower boundary condition of sub-surface hydrology depended on the groundwater level or saturation conditions.The initial conditions of soil water storage were prescribed by setting soil water content.The drainage and deep percolation were calculated based on an empirical linear equation.

Snow dynamics
The snow pack is an important factor to influence the upper boundary condition of soil heat and water flows and has been assumed homogeneous both horizontally and vertically in the model.Precipitation was distinguished between rainfall and snowfall using the temperature as a threshold parameter.Snowmelt was controlled by global radiation, air temperature and heat flux from the soil.The melting caused by global radiation was to some extent dependent on snow age.Liquid water retained in the snow pack was allowed to refreeze.The thermal conductivity of snow was estimated according to snow density.In this study, we used the energy balance approach to calculate snow surface temperature and sensible and latent heat fluxes upon snow surface (Gustafsson et al., 2004).

Plant growth
The potential rate of leaf assimilation was estimated using the light-use efficiency approach, which assumes that the total plant growth is proportional to the absorbed global radiation but limited by unfavourable conditions of Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.temperature, water and nitrogen (Monteith 1972;Monteith and Moss, 1977).Too high or too low leaf temperatures, too high leaf C:N ratio and water stress (i.e. the ratio of transpiration to potential transpiration) could stunt plants' growth and development.The growth stage of plants (e.g.start of growth, leafing and grain development etc.) was determined according to species types and regulated by temperature sums.The allocation of C to different storage compartments (root, stem and leaf) was estimated by the multiplicative response of biomass, the leaf C:N ratio and water stress.The new (the current year) and old (the previous years) C storage compartments were distinguished and had different rates of litterfall.Plant respiration accounting for both maintenance and growth was simulated using the exponential "Q10 type" temperature response function (Hansen and Jensen, 1977).The model represented vegetation structure using the "multiple canopies" approach, which differentiates structure of each canopy stand in their competence of using the common resources, e.g.light, water and nutrients, for their growth.In our case, we parameterized two sets of canopy properties corresponding to the dwarf evergreen shrub (i.e. C. tetragona) and the prostrate forbs (i.e.mosses and lichen) respectively.These two canopies had different sizes in surface cover (C.tetragona: 53%, mosses: 8%), maximum canopy height (C.tetragona: 15 cm, mosses: 5 cm), maximum leaf area index (C.tetragona: 1, mosses: 0.5) and root lowest depth (C.tetragona: 20 cm, mosses: 5 cm).The maximum lifetime of leaves for C. tetragona and mosses were set to 3 years and 1 year respectively.

Soil decomposition
At each soil layer, soil organic matter was stored by a litter pool and a humus pool, which had a high and low C turnover rate correspondingly.The decomposition in both C pools was calculated by following the first order kinetics and constrained by soil temperature and soil moisture response function.The soil temperature response function used "Ratkowsky function", a quadratic function to account for temperature dependency of microbial activities (e.g.bacterial growth) and was controlled by high and low temperature threshold parameters (Ratkowsky, 1982).The soil moisture response function relied on the parameters which control the dependency curve within the high and low interval of soil water content (Skopp et al., 1990).Between these two water content intervals, there is no soil moisture limitation.Since the microbes were implicitly included in the litter pool, the synthesis of microbial biomass and metabolites constituted internal cycling.The initial values for the C and N content in each soil layer were prescribed by measurements and partitioned into the two pools for each layer according to the measured C:N ratio and soil bulk density as described in Table 1.

Model setup, calibration and evaluation
The model used the hourly meteorological measurements (i.e.2m air temperature, precipitation incoming shortwave and longwave radiation, wind speed, and relative humidity) as drivers.By assuming vegetation started to grow from an equilibrium state, we implemented a 90-year run which was repeatedly driven by the meteorological measurements of 1996 as a model spin-up.We applied the GLUE (Generalized Likelihood Uncertainty Estimation, Jansson, 2012)  (III) to assign range and distribution functions for selected parameters; (IV) random realization of the model with independence between parameters and previous runs (Wetterstedt & Ågren, 2011).
In this study, we performed a 19,000 Monte-Carlo multi-run based on the stochastic sampling of 33 parameters (Table 2).The minimum and maximum values of the sampled parameters were given in the Table 2.We used performance indexes (i.e.mean errors -ME and root mean square errors -RMSE) to select around 30 posterior runs based on the calibration against the time-integrated CO2 fluxes (Table 3).It is noted that the acceptance criteria is subjective and attained after we have conducted several initial selection tests.The time-integrated CO2 fluxes refer to the hourly measurements of NEE being transformed to cumulative sequences of fluxes for a day, a year and the entire measurement period.The NEE measurements were not gap-filled.Therefore, the transformations were made to smooth out uncertainties in single measured point and rather emphasize the pattern of flux within the day, within the year or for the long-term trend of the measurements.We calibrated the model based on the transformed NEE fluxes.After calibration, the corresponding behaviour model ensembles were named "the diurnal behaviour model ensemble (DBME), the seasonal behaviour model ensemble (SBME) and the long-term behaviour model ensemble (LBME)".
To identify the importance of controlling parameters and relevant processes reflected by the three behaviour model ensembles, we calculated the Pearson correlation between posterior parameters and R 2 (the coefficient of determination for linear regression) of each measurement variable, and the inter-correlation matrix of parameters.
We use the mean, coefficient of variation (CV, the standard deviation divided by the mean), skewness and kurtosis of probability density function (PDF) and cumulative distribution of frequency (CDF) to describe how well the posterior parameters have been constrained by each behaviour model ensemble.Finally, the mean and trend for the inter-annual variation of internal C fluxes and C storage pools in the three behaviour model ensembles were presented.

Spectral analysis
In order to investigate how well the three behaviour model ensembles were constrained to reflect multiple time scales of ecosystem behaviours, we calculated a multiscale variance of "maximal overlap discrete wavelet transform" (MODWT) for the measured and simulated CO2 fluxes.This approach decomposes the data sets into time-varying components to reveal the data patterns within single time window.The wavelet analysis was conducted using the Matlab wavelet toolbox on the basis of the "Haar" mother wavelet, which has been widely used to analyse the atmospheric turbulence data (Hudgins et al., 1993).We firstly pre-processed the data sets by replacing all the missing values in the measurements with 0 so as to have a complete time series.The model results only selected the values within the measurement dates and used 0 for the days with no measurements.The wavelet decomposition was performed for 15 levels.Each level corresponds to different time window length with the power of 2 at hours.

Wavelet detection of time-varying variance for the measurements and model ensemble means
The variance of MODWT for the measured and modelled CO2 fluxes was presented over the time scales in power of 2 for hours based on 15 transformed levels (Fig. 2).Accordingly, the patterns of the variance changing from the diurnal scale (hours) to the inter-annual scale (a year, equivalent to 2 13.09 hours) were well distinguished for each time series of dataset.The behaviour model ensemble means generally agreed well with the measurements in the patterns of variance across a wide range of time scales, but differed a bit in the magnitude at some individual time windows.Within the daily to weekly time scale (i.e. from 2 2 to 2 7 hours), DBME showed fewer deviations from the measurements than the other two model ensembles, and LBME showed the largest deviation from the measurements.On the contrary, at a relatively long time scale (e.g. 2 9 to 2 11 hours), the least deviation between the measurements and the model was found in LBME.

The model performance in simulating the measurement variables
The posterior runs effectively constrained the prior NEE fluxes by narrowing down the range of ME and RMSE (Fig. 3).All the behaviour model ensembles displayed a symmetrical Gaussian distribution of ME in NEE centring on 0 (Fig. 3a).They also showed their capabilities to best describe the time-scale integrated C fluxes based on their defined behaviours (Fig. 3b-3c).For instance, the smallest RMSE for the daily-integrated NEE was found in DBME (Fig. 3b), and this is similar to SBME and LBME, which had the smallest RMSE in the yearly-integrated NEE and the entire period-integrated NEE respectively (Fig. 3c and 3d).Despite that the model used the observed NEE as the only constraint, the other measurement variables were fairly well simulated (Fig. 4).It implies the holistic ecosystem dynamics represented by the model is reasonably consistent.In each ensemble, R 2 of most measurement variables was larger than 0.5 and the averaged ME was close to 0. However, the model showed a relatively large underestimation on soil water content (SW0.1m and SW0.3m) and latent heat flux (LHF) (Fig. 4a).The relatively low R 2 in NEE in DBME was partly due to the large uncertainties in simulating the autumn-time NEE (NEEOct), which had a lower R 2 than NEE and NEE of growing season (NEEJul-Oct).It is noted that R 2 of the calibrated variables is sensitive to the length of a cumulative time sequence so that R 2 became larger as the cumulative time sequence was lengthened.Although all the ensembles had overestimated snow depth, they still simulated reasonable soil temperatures by using the posterior parameters.

The model residuals allocated at the diurnal, seasonal and long-term time scales
The residuals between the model and measurements were allocated for the diurnal, seasonal and long-term course (Fig. 5).The diurnal patterns of the hourly-mean residuals showed that the model overestimated C uptake during the daytime and C release during the night-time.LBME showed a wider uncertainty band and greater amplitude of the residuals than the other two ensembles (Fig. 5a).For the seasonal patterns of the daily mean residuals, the model overestimated C uptake in early spring and later autumn and C release in summer, and DBME and SBME had a Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.lower amplitude of deviations than LBME (Fig. 5b).In autumn, SBME showed fewer residuals than DBME.For the total cumulative residuals, LBME was constrained much better than the other two ensembles (Fig. 5c).In particular, DBME seemed to have divergent trends for the inter-annual variation of the cumulative residuals.For the period 2000-2007, a positive bias was propagated, but for the period 2008-2014, a negative bias tended to offset the positive bias and resulted in a good balance of the two at the end of the period.

Inter-annual and seasonal variation of NEE
All the model ensembles showed a similar inter-annual and seasonal variation of NEE, but they also differed in the slope of the linear trend and the year-round C budget (Fig. 6).The annual mean NEE started from a positive value (C source) at 2000 and varied towards a negative value (C sink) in most years (Fig. 6a).All of the ensembles reached the highest C uptake in 2008.After 2008, there was no trend for DBME, but both SBME and LBME indicated that the annual NEE tended to shift towards a C source.However, for the entire period, only DBME showed a significant downward trend of NEE (Table 4).For the seasonal profile of NEE, the NEE started from a C source in early spring and became an increased C sink in the growing season, and then was offset by increased C release in autumn (Fig. 6b).DBME showed the highest C uptake in summer and the lowest C release in spring and autumn.The year-round C budget for the three ensembles showed the ecosystem approximated to either a week sink (i.e.DBME: -9.8 ± 4.4 g C m -2 yr -1 and SBME: -4.1 ± 4.7 g C m -2 yr -1 ) or a neutral C balance (LBME: -0.1 ± 6.4 g C m -2 yr -1 ).

The mean C budget and the temporal trends for C fluxes and stores based on the behaviour model ensemble mean
In general, compared to the long-term behaviour model ensemble, the short-term behaviour model ensemble showed a lower estimate for the mean C flux rates and the plant-associated C storage pools and a higher estimate for the soil C pools (Table4).All the three behaviour model ensembles displayed a slightly increasing trend for the ecosystem C uptake (i.e. the negative NEE trend), but only the trend of DBME was significant (Table 4).Other significant trends in DBME were found in the inter-annual variation of GPP, ER, plant respiration (Rplant), root respiration (Rroot), litterfall, and all the C pools but the litter pool (Clitter).By contrast, LBME had the fewest number of C fluxes with significant trends.The GPP and ER simulated by DBME were -59.4 ± 16.1 g C m -2 yr -1 and 41.3 ± 4.4 g C m -2 yr -1 respectively, around 60% of the C flux rates estimated by LBME.The difference of C balance in the three model ensembles was largely subject to the proportion of soil CO2 production.The fractions of respiration for humus and soil (i.e.Rhumus and Rsoil) in DBME, SBME and LBME were 20.8%, 27.8% and 33.3% respectively.A higher portion of soil respiration simulated by LBME originated from changes in the soil humus pool (Chumus), indicating a larger amount of old C release.Opposed to the trend of Clitter, the trend of Chumus was found significant for all the ensembles, with a higher magnitude for a longer term behaviour model ensemble.It is noted that the moss C pool (Cmoss) was higher in DBME than that in SBME and LBME.

Sensitivities of posterior parameters to model performance
We counted the number of posterior parameters with the Pearson correlation coefficient for linear regression (r > 0.3 or r < -0.3) to the posterior R 2 of the measurement variables and grouped them into the processes they belonged to.This helped us to identify which processes were important to which measurement variables and how this indication differed among the model ensembles.All the model ensembles showed a certain pattern on the count of sensitive parameters for their posterior values and R 2 of the measurement variables (Fig. 7).The parameters were grouped into the processes of snow dynamics, soil water and soil heat, plant growth, and soil decomposition.The surface temperature (ST0m) was highly associated with the parameters in the processes of snow dynamics, soil heat and soil water and plant growth.But for soil temperature at the depth of 0.6 m and 1.3 m (i.e.ST0.6m and ST1.3m), LBME showed the sensitive parameters belonged to the processes of soil decomposition.For the soil water content at 0.1 m (SW0.1m),parameters from all the processes were sensitive to the model performance, but this differed from the soil water content at 0.3 m (SW0.3m), which may exclude parameters from one or another process.The sensitive parameters to snow depth common to all the model ensembles were from the processes of plant growth and snow dynamics.For the radiation fluxes, the common sensitive parameters were from the processes of plant growth, snow dynamics and soil water and soil heat.The parameters in the process of soil decomposition appeared sensitive to latent heat flux in both SBME and LBME.For NEE of the entire period or sub-periods, the soil decomposition parameters increased their number in the longer time behaviour models.Overall, the performance of SBME for most measurement variables was correlated to the parameters from all the four processes, while in DBME, more correlations were found in the parameters associated with snow dynamics and soil water and heat transfer, and in LBME, more correlations were found in the parameters associated with soil decomposition.Most parameters associated with plant growth showed higher sensitivities than other parameters in all the three model ensembles.

Inter-correlations between posterior parameters
The inter-correlations between posterior parameters reflected how close the inter-links of processes were controlled by the parameters in the behaviour models.This may infer that the data used to constrain the model was not suitable to distinguish various possible explanation to the measurements or that the ecosystem had an inherent dependence between different properties.The Pearson correlation coefficient for linear regression (r > 0.3 or r < -0.3) was adopted to indicate a moderate or high correlation between parameters.Among all the model ensembles, DBME showed that parameters with a higher number of significant correlations to others were from the processes of snow dynamics; plant respiration and soil decomposition (Fig. S3a, e.g.DensityCoefWater (4), AlbSnowMin (4), EquiAdjustPsi (4), RootRate (4), RateCoefLitter (4), RateCoefHumus (5), TempMin (6); it is noted that the number within the brackets indicates the number of parameters with a significant correlation).SBME showed that parameters with a higher number of significant correlations to others were relevant to snow density, surface evaporation, surface temperature, and plant growth (Fig. S3b, e.g.DensityOfNewSnow (4), EquiAdjustPsi (5), CFrozenSurfCorr (5), TLMin (4), rOptimum_value (6), FixNsupply (6), MCoefLeaf (5)).LBME revealed more number of significant inter-correlations for parameters were related to the processes of soil heat and water and plant growth (Fig. S3c, e.g.RoughLMomSnow (5), DensityCoefMass (5), EquiAdjustPsi (4), FreezepointFWi (6), TLMin (7), TLOpt1(5), TempMin (5)).In general, we conclude that DBME showed more unconstrained parameters in soil decomposition, which however showed the best constraint in LBME.SBME revealed more inter-correlation for parameters in plant growth, indicating that more equifinalities may exist in this posterior ensemble.

Posterior distribution of calibrated parameters
Most parameters showed a negative kurtosis indicating that the posterior PDF of parameters had a larger combined weight of the tails relative to the rest of the distribution, and thus their distribution had a flatter peakedness comparing to the normal distribution (Table 2).Few parameters (e.g.RoughMomSnow in DBME and LBME, EquilAdjustPsi in DBME, TLOpt1 in LBME, TLMax in SBME and RateCoefHumus in DBME) showed a relatively high skewed distribution (i.e.skewness is <-1 or >1).In the process of snow dynamics, RoughLMomsnow and AlbSnowMin were better constrained than other parameters (Fig. 8).In the process of soil heat and water transfer, all the model ensembles provided a similar posterior PDF for the parameter OrganicLayerThick and larger differences for the parameters EquilAdjustPsi, CFrozenSurfCorr, SurfCoef and FreezepointFWi.In the process of plant growth, a better constraint on the posterior PDF in the model ensembles was found on the parameters related to the temperature response function of photosynthesis (i.e.TLMin, TLOpt1, TLOpt2 and TLMax) and plant respiration (MCoefLeaf, MCoefRoot, ResptemQ10 and rOptimum).The parameters governing the litterfall rates of leaf and root (i.e.LeafRate and RootRate) showed a similar constraint in all the three model ensembles.In the process of soil decomposition, the model ensembles displayed distinguished patterns for all the parameters, for instance, the rate coefficient of humus (RateCoefHumus) suggested much lower and higher posterior means by DBME and LBME respectively.subtracting the daytime respiration.Their growing season GPP and ER for 2000-2010 were -78.6 ± 13.7 g C m -2 yr -1 and -57.9 ± 11.4 g C m -2 yr -1 respectively.In contrast, this study is based on a process-oriented model to estimate the yearly budget of NEE.For other sites in pan-Arctic, the year-round measurements in Svalbard and Alaska all indicate that the current tundra ecosystem is now in the transition from a C sink to a C source (e.g.Oechel et al., 2014;Lüers et al., 2014 andEuskirchen et al., 2016).Particularly, Ecskirchen et al. ( 2016) reported a C. tetragona heath ecosystem in the Alaskan Arctic tundra had a yearly NEE around 20 g C m -2 for 2014-2015, which was largely affected by the freezing time of active layers.In our study, LBME showed a larger tendency for ecosystems to become a C source by suggesting a higher amplitude for the diurnal CO2 flux and a higher rate of humus decomposition.This implies that a faster C cycling, a higher rate of CO2 release associated with the old C decomposition rate and the thawing depth can be possible drivers for the ecosystem sink-source transition.

Model performance and parameter uncertainties
The model-measurement residuals appearing at the same time scale across all the model ensembles imply that the model may need more measurement variables as additional constraints or examine the model structure or parameterization.For instance, to address the bias in the diurnal cycle, the sub-daily soil temperatures should be included in the calibration.For the diurnal NEE, we found a negative bias occurred at noon and a positive bias occurred in the evening.Both biases are likely associated with the absorbed light in the early spring and midsummer (Fig. S1).In the model, the simulated leaf area and snow coverage are important factors to estimate the absorbed light correctly.In the seasonal cycle, the negative bias in the early growing seasons can be attributed to the earlier photosynthesis occurring in the model (Fig. S2).Some studies (e.g.Mä kelä et al., 2004;Mellander et al., 2008;Wu et al., 2012b) showed that boreal ecosystems may encounter the effects of thermal acclimation for photosynthesis in the early spring, which can delay the start of the photosynthesis.Soil moistures should be an additional constrain for the growing season GPP and ER, because the freezing and thawing of active layers is highly dependent on the soil moisture.To address the bias of NEE occurring in the long-term time scale, the biophysical properties of plants should be included as the model calibration variables or the model drivers.
The behaviour model ensembles showed different emphasis on the specific processes which affect the performance of the model to simulate measurement variables, but the plant growth processes were the common processes having more sensitive parameters than other processes in all the ensemble runs.The posterior uncertainties of many parameters were substantially reduced, and some of them show the distinguished patterns in the three model ensembles.The exponential coefficient in the Q10 equation for plant respiration were suggested a higher value range (2.75-3.5)than the default value (0.5 -1.75).The default value was often applied to the boreal forests, which had warmer ambient temperature.Higher Q10 indicates the temperature sensitivities is lower in the lower temperature but higher in the warm temperature (i.e.20 °C) in our study.All the model ensembles suggested a lower temperature as the cold threshold for temperature response of soil respiration, indicating a convergence of temperature sensitivity.
Further, there may have some equifinalities for the posterior parameters, which allow different parameter sets to provide equal efficiency to describe the system, and suggest right estimates because of wrong reasons.This issue can be addressed by adding more independent data to constrain the model (William et al., 2009;Carvalhais et al., 2010).
Page 12 of 34

Important drivers for the diurnal, seasonal and long-term behaviours of ecosystem CO2 exchange in high Arctic tundra
To elucidate important divers for the temporal behaviours of ecosystem CO2 exchange, we presented the diagnostic variables (i.e. the response anomaly) to indicate how photosynthesis and ER are regulated by their responses to air temperature and water stress, absorbed radiation and responses of litter and humus decomposition varying at the daily, seasonal and long-term time scales.The responses of litter and humus decomposition have accounted for depth-specific soil temperature, soil moisture and substrate concentration.The largest diurnal deviation in the absorbed radiation was found highest at noon but lowest at night (Fig. 9).The temperature response showed positive anomaly in the daytime with a larger portion in the afternoon.The response of humus and litter respiration occurred mostly in the afternoon, whereas soil moisture response was much lower in the daytime than night-time.For the seasonal profile, the highest temperature response occurred in the mid-growing season, and the absorbed radiation started earlier in conjunction with the snowmelt period.The humus and litter respiration also had a higher responses in the mid-summer but humus respiration response seemed to have a little delayed than litter response, possibly because of the delayed response of deeper soil horizons where old carbon has a larger proportion.The soil moisture response decreased in the growing season, and a larger decrease was seen in autumn.For the long-term scale, the response of litter respiration and soil moisture seem to level off for the entire period, while the humus respiration response had a little reduction in the last few years.The year-to-year variation of litter and humus respiration response is largely in line with the tendency of temperature response.The absorbed radiation and temperature response showed a large increase before 2008 and level off for the rest the period.The absorbed radiation to some extent reflected the biological properties response.For instance, the surface albedo has a high decrease in the first half period than the second half period because of changes in vegetation cover.Note that all the above-mentioned responses were quantified by the corresponding time scale behaviour model ensemble means.

Implication for both modelling and measurements
Our study demonstrates a clear trade-off when the model is tuned to capture both the short-and long-term patterns in ecosystem CO2 exchange.This agrees the previous findings that environmental and biotic factors might represent different roles in explaining fluctuations in CO2 fluxes across time scales (Richardson et al. 2007;Wu et al., 2012a;Wu et al., 2017).However, as we have not validated the model with any observation that reflects seasonal and interannual variabilities of biological responses, for instance, canopy cover or LAI, we cannot conclude that the large bias in the total cumulative CO2 fluxes in DBME is attributed to uncertainties in model parameterization or model structure deficiency.The former uncertainties can be addressed by using observations to prescribe the dynamics of parameters, (e.g.LAI).The latter uncertainties can be addressed by including processes which describe a more appropriate seasonal pattern of ecosystem dynamics.For instance, the overestimated C uptake appeared frequently in spring for several years.This implies that the model may have to account for thermal acclimation effects of photosynthesis to explain the delay of C uptake.LBME suggests a more precise description of long-term behavior of CO2 fluxes by estimating a faster turnover rate of soil C pool.This also highlights that to quantify the turnover rate of soil C is crucial to project the long-term trend of ecosystem CO2 exchange.
The high frequency of measurements used in this study allows the process-oriented model to identify the main drivers for variation of C fluxes across the wide range of time scales and to quantify the year-around C budget.The process-oriented modelling with a strong emphasis on the balance of energy, water and carbon at a high resolution for soil profile and temporal evolution has demonstrated the efficiency in process representation and the challenges of capturing both short-term and long-term CO2 variation.The non-growing season respiration and the burst events of CO2 efflux in spring in relation to seasonal and future climate trends calls for the extension of monitoring campaign to cover the entire non-growing-seasons.Moreover, it is important to investigate to which extent the early spring burst originates from the winter CO2 production and the partition of soil CO2 production from the decomposition of labile C storage and old C storage located at the near-surface horizon.To compromise the uncertainties of winter precipitation, snow properties like snow cover extent and snow water equivalent are of high interest in that they determine soil temperature and moisture conditions for the non-growing season respiration to regulate the current transition of NEE.Further studies are needed to investigate plant growth (i.e.photosynthesis and respiration) in the subnivean microclimate.

Conclusions
This study demonstrates that three behaviour model ensembles constrained by the time-integrated CO2 fluxes using the Monte-Carlo runs were able to describe the variation of 15 years of eddy covariance NEE measurements for a high Arctic heath ecosystem on a daily, seasonal and long-term basis.The inter-annual variation of NEE showed a trend from the ecosystem being a C sink to a C neutral balance for all the three behaviour model ensembles.The long-term behaviour model ensemble simulated a more intensified diurnal C cycle than the short-term behaviour model ensemble.The intensified C cycle was mainly attributed to a faster depletion of soil C pools and higher amplitude of diurnal CO2 variation.The correlations of posterior parameters and R 2 reflected that parameters in the processes of soil water and heat and snow dynamic regulated the short-term behaviour of CO2 fluxes, while parameters in soil decomposition processes regulated the long-term behaviour of CO2 fluxes.Our results suggest that using the time-integrated CO2 fluxes as model constraints can be a good diagnostic approach to evaluate how the performance of the model is appropriate for different time scales of processes and if there is any transition in ecosystem processes the model misrepresents.However, more efforts in quantifying models' uncertainties and more independent measurements are still needed to further improve our understanding the key drivers of high Arctic ecosystems C dynamics associated with substantial changes observed in the environmental conditions.Moss biomass (g m -2 ) 235 a or 120 b g m -2 Leaf nitrogen (g m -2 ) 1 g m framework to calibrate the model.The GLUE framework adopts a Monte-Carlo based approach to reduce parameter uncertainties of the prior runs and thus attain the behaviour model ensembles of accepted runs.It normally involves the following procedures: (I) to define the acceptance criteria; (II) to select uncertain parameters; Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.
Our model ensembles demonstrate the strength of describing the variability of C fluxes at diurnal, seasonal and long-term time scales.The wavelet analysis further supports that the model constrained using the time-integrated CO2 flux has effectively distinguished the model behaviour by matching the observational variance across discrete time windows.The multi-yearly mean of NEE for 2000-2014 represented by the model ensemble means ranges from -9.8 ± 10.4 g C m -2 of DBME to -0.1 ± 10.5 g C m -2 yr -1 of LBME, which are higher than the growing season NEE measurements reported by previous studies for the same site (e.g.-12.4 ± 8.7 g C m -2 yr -1 for the period of 1997 and 2000-2003 by Grøndahl et al. (2007) and -20.6 ± 11 g C m -2 yr -1 for 2000-2014 by Lund et al. (2012)).The large year-to-year variability of NEE cannot distinguish the ecosystem as a significant C sink or source, and the uncertainties of ensemble runs can even further enlarge the year-to-year variability of NEE.Many other modelling studies (e.g.McGuire et al. (2012); Fisher et al. (2014); Zhang et al. (2013a)) also showed that the process-based models have simulated the present-day tundra as a week C sink, but most of these models have only been calibrated by the growing season NEE measurements.Lund et al. (2012) used a light response curve model to estimate GPP by Page 11 of 34 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .Figure 2 .
Figure 1.Physical/physiological controls and biological responses operated at different temporal (x axis) and ecosystem-

Figure 3 .
Figure 3. (a) The number of runs in mean errors (ME) of NEE flux (g C m -2 ).(b) Cumulative distribution of frequency

Figure
Figure 4.The coefficient of determination for linear regression (R 2 ) and mean errors (MR) for the behaviour model

Figure 8 .
Figure 8. Cumulative distribution frequency (CDF) of all the sampled parameters in the prior runs (black) and posterior

Figure 9 .
Figure 9.The normalized anomalies of ecosystem responses to total photosynthesis and ecosystem respiration at the minimum and maximum parameter values in the prior runs.The mean (Mean), coefficient of variation (CV, the standard deviation divided by the mean), , https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.Table 3.The acceptance criteria used to constrain a 19,000 Monte-Carlo multi-run for three behaviour model ensembles.mean C budget simulated by the three model behaviour ensembles (DBMEdiurnal behaviour model ensemble, SBMEseasonal behaviour model ensemble and LBMElong-term behaviour model ensemble).The bold type indicates the slope of linear trend with p Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-382Manuscript under review for journal Biogeosciences Discussion started: 29 September 2017 c Author(s) 2017.CC BY 4.0 License.

Table 2 .
The