Evaluation of terrestrial pan-Arctic carbon cycling using a data-assimilation system

There is a significant knowledge gap in the current state of the terrestrial carbon (C) budget. Recent studies have highlighted a poor understanding particularly of C pool transit times and of whether productivity or biomass dominate these biases. The Arctic, accounting for approximately 50 % of the global soil organic C stocks, has an important role in the global C cycle. Here, we use the CARbon DAta MOdel (CARDAMOM) data-assimilation system to produce pan-Arctic terrestrial C cycle analyses for 2000–2015. This approach avoids using traditional plant functional type or steady-state assumptions. We integrate a range of data (soil organic C, leaf area index, biomass, and climate) to determine the most likely state of the high-latitude C cycle at a 1^∘ × 1^∘ resolution and also to provide general guidance about the controlling biases in transit times. On average, CARDAMOM estimates regional mean rates of photosynthesis of 565 g C m⁻² yr⁻¹ (90 % confidence interval between the 5th and 95th percentiles: 428, 741), autotrophic respiration of 270 g C m⁻² yr⁻¹ (182, 397) and heterotrophic respiration of 219 g C m⁻² yr⁻¹ (31, 1458), suggesting a pan-Arctic sink of −67 (−287, 1160) g Cm⁻² yr⁻¹, weaker in tundra and stronger in taiga. However, our confidence intervals remain large (and so the region could be a source of C), reflecting uncertainty assigned to the regional data products. We show a clear spatial and temporal agreement between CARDAMOM analyses and different sources of assimilated and independent data at both pan-Arctic and local scales but also identify consistent biases between CARDAMOM and validation data. The assimilation process requires clearer error quantification for leaf area index (LAI) and biomass products to resolve these biases. Mapping of vegetation C stocks and change over time and soil C ages linked to soil C stocks is required for better analytical constraint. Comparing CARDAMOM analyses to global vegetation models (GVMs) for the same period, we conclude that transit times of vegetation C are inconsistently simulated in GVMs due to a combination of uncertainties from productivity and biomass calculations. Our findings highlight that GVMs need to focus on constraining both current vegetation C stocks and net primary production to improve a process-based understanding of C cycle dynamics in the Arctic.

1 Introduction

Arctic ecosystems play a significant role in the global carbon (C) cycle (Hobbie et al., 2000; McGuire et al., 2012). Slow organic matter decomposition rates due to cold and poorly drained soils in combination with cryogenic soil processes have led to an accumulation of large stocks of C stored in the soils, much of which is currently held in permafrost (Tarnocai et al., 2009). The permafrost region soil organic C (SOC) stock is more than twice the size of the atmospheric C stock and accounts for approximately half of the global SOC stock (Hugelius et al., 2014; Jackson et al., 2017). High-latitude ecosystems are experiencing a temperature increase that is nearly twice the global average (AMAP, 2017). The expected future increase in temperature (IPCC, 2013) and precipitation (Bintanja and Andry, 2017) will likely have consequences for the Arctic net C balance. As high latitudes warm, C cycle dynamics may lead to an increase in carbon dioxide (CO₂) emissions through ecosystem respiration (R_eco) driven by, for example, larger heterotrophic respiration (Commane et al., 2017; Schuur et al., 2015; Zona et al., 2016), drought stress on plant productivity (Goetz et al., 2005), and episodic disturbances (Lund et al., 2017; Mack et al., 2011). Alternatively, temperature-induced vegetation changes (Forkel et al., 2016; Graven et al., 2013; Lucht et al., 2002) may increase gross primary productivity (GPP) due to extended growing seasons (Zeng et al., 2011), CO₂ fertilisation (Zhuang et al., 2006), and shifts in vegetation cover such as greening (Myneni et al., 1997; Zhu et al., 2016) and shrub expansion (Myers-Smith et al., 2011). Consequently, ecosystem responses may feed back on climate with unclear magnitude and sign (Anav et al., 2013; Murray-Tortarolo et al., 2013; Peñuelas et al., 2009). As a result of the significant changes that are already affecting the structure and function of Arctic ecosystems, it is critical to understand and quantify the historical C dynamics of the terrestrial tundra and taiga and their sensitivity to climate (McGuire et al., 2012).

Although the land surface is estimated to offset ∼30 % of anthropogenic emissions of CO₂ (Canadell et al., 2007; Le Quéré et al., 2018), the terrestrial C cycle is currently the least constrained component of the global C budget and large uncertainties remain (Bloom et al., 2016). Despite the importance of Arctic tundra and taiga biomes in the global land C cycle, our understanding of interactions between the allocation of C from net primary productivity (NPP), C stocks (C_stock), and transit times (TTs) is deficient (Carvalhais et al., 2014; Friend et al., 2014; Hobbie et al., 2000). The TT is a concept that represents the time it takes for a particle of C to persist in a specific C stock, and it is defined by the C stock and its outgoing flux, here addressed as TT = C_stock∕NPP. According to a recent study by Sierra et al. (2017), TT is an important diagnostic metric of the C cycle and a concept that is independent of model-internal structure and theoretical assumptions for its calculation. Terms such as residence time (Bloom et al., 2016; Friend et al., 2014), turnover time (Carvalhais et al., 2014), and turnover rate (Thurner et al., 2016; TT = 1∕turnover rate) are used in the literature to represent the concept of TT (Sierra et al., 2017). Studies have focused more on the spatial variability with climate of ecosystem productivity rather than C transit times (Friend et al., 2014; Nishina et al., 2015; Thurner et al., 2016, 2017). Friend et al. (2014) detailed that transit time dominates uncertainty in terrestrial vegetation responses to future climate and atmospheric CO₂. They found a 30 % larger variation in modelled vegetation C change than response of NPP. Nishina et al. (2015) also suggested that long-term C dynamics within ecosystems (vegetation turnover and soil decomposition) are more critical factors than photosynthetic processes (i.e. GPP or NPP). The respective contribution of bias from biomass and NPP to biases in transit times remains unquantified. Without an appropriate understanding of the current state and dynamics of the C cycle, its feedbacks to climate change remains highly uncertain (Hobbie et al., 2000; Koven et al., 2015).

There are currently efforts to incorporate both in situ and satellite-based datasets to assess C cycle retrievals and to reduce their uncertainties. At a local scale, the net ecosystem exchange (NEE) of CO₂ between the land surface and the atmosphere is usually measured using eddy covariance (EC) techniques (Baldocchi, 2003). International efforts have led to the creation of global networks such as FLUXNET (http://fluxnet.fluxdata.org/, last access: 9 April 2019) and ICOS (https://www.icos-ri.eu/, last access: 9 April 2019), to harmonise data and support the reduction of uncertainties around the C cycle and its driving mechanisms. However, upscaling field observations to estimate the regional to global C budget presents important challenges due to insufficient spatial coverage of measurements and heterogeneous landscape mosaics (McGuire et al., 2012). Furthermore, harsh environmental conditions in high-latitude ecosystems and their remoteness complicates the collection of high-quality data (Lafleur et al., 2012). Given the lack of continuous, spatially distributed in situ observations of NEE in the Arctic, it remains a challenging task to calculate with certainty whether or not the Arctic is a net C sink or a net C source and how the net C balance will evolve in the future (Fisher et al., 2014). Over the past decade, regional to global products generated from in situ networks and/or satellite observations have improved our understanding of the terrestrial C dynamics. These range from machine-learning-based upscaling of FLUXNET data (Jung et al., 2017), remotely sensed biomass products (Carvalhais et al., 2014; Thurner et al., 2014), and the creation of a global soil database (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012). Due to a reliance on interpolation and upscaling with other spatial data, it is challenging to evaluate these products for inherent biases.

Global vegetation models (GVMs) have been developed to determine global terrestrial C cycling, through representing vegetation and soil processes, including vegetation dynamics (i.e. growth, competition, and turnover) and biogeochemical (i.e. water, carbon, and nutrients cycling) responses to climate variability (Koven et al., 2011; Sitch et al., 2003; Woodward et al., 1995). The advantage of using process-based models to characterise C dynamics is that processes which drive ecosystem–atmosphere interactions can be simulated and reconstructed when data are scarce. However, C cycle modelling in GVMs typically relies on parameters retrieved from literature, prescribed plant functional type (PFT), and a spin-up process ensuring C stocks (biomass and SOC) reach steady state. Further, inherent differences of model structure contribute more significantly to GVM uncertainties (Exbrayat et al., 2018; Nishina et al., 2014) than do differences in climate projections (Ahlström et al., 2012). Many model inter-comparison projects have demonstrated a lack of coherence in future projections of terrestrial C cycling (Ahlström et al., 2012; Friedlingstein et al., 2014). Recent studies have used simulations from the first phase of the Inter-Sectoral Impact Model Inter-comparison Project (ISI-MIP) (Warszawski et al., 2014) to evaluate the importance of key elements regulating vegetation C dynamics but also the estimated magnitude of their associated uncertainties (Exbrayat et al., 2018; Friend et al., 2014; Nishina et al., 2015; Thurner et al., 2017). An important insight is that TTs in GVMs are a key uncertain feature of the global C cycle simulation. Further, GVMs tend not to report uncertainties in their estimates of stocks and fluxes, which weakens their analytical value.

To address these issues we integrate model and data more formally. We apply data assimilation (DA), defined as a Bayesian calibration process for a model of a dynamic system. DA, through probabilistic parameterisation, supports robust model estimates of C stocks and fluxes consistent with multiple observations and their errors (Fox et al., 2009; Luo et al., 2009; Williams et al., 2005). By following Bayesian methods, the uncertainty in observations weights the degree of data constraint, and the outcome is a set of acceptable parameterisations for a given model structure linked to likelihoods. Overall, this approach determines whether model structure, observations, and forcing are (in)consistent and thus assesses the validity of model structure. By assimilating co-located climatic, ecological, and biogeochemical data from remote-sensing observations at a specific grid scale across landscapes and regions, DA can map parameter estimation and uncertainties.

Here, we use the CARbon DAta MOdel framework (CARDAMOM) (Bloom et al., 2016; Bloom and Williams, 2015; Smallman et al., 2017) to analyse the pan-Arctic terrestrial carbon cycle at 1^∘ resolution for the 2000–2015 period. We assimilate gridded observations of leaf area index (LAI), biomass, and SOC stocks at these spatio-temporal scales into an intermediate-complexity C model (DALEC2, which is less complex than GVMs). We compare analyses of C dynamics of Arctic tundra and taiga with (a) global products of GPP (Jung et al., 2017) and heterotrophic respiration (R_h) (Hashimoto et al., 2015); (b) NEE, GPP, and R_eco field observations from eight high-latitude sites included in the FLUXNET2015 dataset; and (c) six GVMs from the ISI-MIP2a comparison project (Akihiko et al., 2017). Our objectives are to (1) present and evaluate the analyses and uncertainties of the current state of the pan-Arctic terrestrial C cycling using a DA system, (2) quantify the degree of agreement between the CARDAMOM product with local- to global-scale sources of available data to assess analytical bias, and (3) use CARDAMOM as a benchmarking tool for the ISI-MIP2a models to provide general guidance towards GVM improvements in transit time simulation. Finally, we suggest future work to be done in the context of advancing pan-Arctic C cycle modelling.

2 Data and methods

2.1 Pan-Arctic region

The spatial domain we considered in this study (Fig. S1 in the Supplement) corresponds to the extent of the Northern Circumpolar Soil Carbon Database version 2 (NCSCDv2) dataset (Hugelius et al., 2013a, b), bounded by 42–80^∘ N and 180^∘ W–180^∘ E and at a spatial resolution of 1^∘ × 1^∘. This area of study totals 18 000 000 km² of land area. We used the GlobCover vegetation map product developed by the European Space Agency (Bontemps et al., 2011) to separate regions dominated by non-forested (tundra) and forested (taiga) land cover types. A complete description of the classes included in each domain can be found in Fig. S1 and caption. The differentiation between tundra and taiga grid cells is in agreement with the tree line delimitated by Brown et al. (1997) together with the tundra domain defined from the Regional Carbon Cycle Assessment and Processes Activity reported by McGuire et al. (2012). The extensive grasslands without presence of trees in some areas such as in south Russia, Mongolia, and Kazakhstan were neglected to focus on higher latitudes. This classification of tundra and taiga totals 8 100 000 and 9 900 000 km² of land area, respectively.

2.2 The CARbon DAta MOdel framework

Here we use CARDAMOM (Bloom et al., 2016) (list of acronyms can be found in Table S1 in the Supplement) to retrieve terrestrial C cycle dynamics, including explicit confidence intervals, in the pan-Arctic region. CARDAMOM consist of two key components: (1) an ecosystem model, the Data Assimilation Linked Ecosystem Carbon version 2 (DALEC2) (Bloom and Williams, 2015; Williams et al., 2005), constrained by observations, and (2) a data-assimilation system (Bloom et al., 2016). This framework reconciles observational datasets as part of a representation of the terrestrial C cycle in agreement with ecological theory.

2.2.1 DALEC2

The DALEC2 ecosystem model simulates monthly land–atmosphere C fluxes and the evolution of six C stocks (foliage, labile, wood, roots, soil organic matter – SOM – and surface litter) and corresponding fluxes. DALEC2 includes 17 parameters controlling the processes of plant phenology, photosynthesis (GPP), allocation of primary production to respiration and vegetation carbon stocks and plant and organic matter turnover rates, all established within specific prior ranges based on ecologically viable limits (Table S2; most priors are uniform with broad ranges). DALEC2 simulates canopy-level GPP via the Aggregated Canopy Model (ACM; Williams et al., 1997), and the most sensitive ACM parameter, related to canopy photosynthetic efficiency, is included in the CARDAMOM calibration. DALEC allocates net primary production to the four plant stocks (foliage, labile, wood, and roots) and autotrophic respiration (R_a) as time-invariant fractions of GPP. Plant C decays into litter and soil stocks where microbial decomposition generates heterotrophic respiration (R_h). Turnover of litter and soil stocks is simulated using temperature-dependent first-order kinetics. For practical purposes we aggregated the different C stocks into photosynthetic (C_photo; leaf and labile), vegetation (C_veg; leaf, labile, wood, and roots), soil (C_dom; litter and SOM), and total ($C_{\mathrm{tot}}=C_{\mathrm{photo}}+C_{\mathrm{veg}}+C_{\mathrm{dom}}$) C stocks. The NEE is calculated as the difference between GPP and the sum of the respiration fluxes ($R_{\mathrm{eco}}=R_{\mathrm{a}}+R_{\mathrm{h}}$), while NPP is the difference between GPP and R_a. Only NEE follows the standard micrometeorological sign convection presenting the uptake of C as negative (sink) and the release of C as positive (source); both GPP and R_eco are reported as positive fluxes. In this study, we addressed C turnover rates and decomposition processes as their inverse rates; this is the C transit time (TT_photo, TT_veg and TT_dom), represented as the ratio between the mean C stock and the mean C input into that stock during the simulation period.

2.2.2 Data-assimilation system

The intermediate complexity of the DALEC2 model compared to typical GVMs facilitates computationally intense Monte Carlo (MC) data-assimilation to optimise the initial stock conditions and the 17 process parameters that shape C dynamics. CARDAMOM is forced with climate data from the European Centre for Medium-Range Weather Forecast Reanalysis interim (ERA-interim) dataset (Dee et al., 2011) monthly for the 2000–2015 period. A Bayesian Metropolis-Hastings Markov chain MC (MHMCMC) algorithm is used to retrieve the posterior distributions of the process parameters according to observational constraints and ecological and dynamic constraints (EDCs; Bloom and Williams, 2015). EDCs ensure that DALEC2 simulations of the terrestrial carbon cycle are realistic and ecologically viable and help to reduce the uncertainty in the model parameters by rejecting estimations that do not satisfy different conditions applied to C allocation and turnover rates as well as trajectories of C stocks.

Observational constraints include monthly time series of LAI from the MOD15A2 product (Myneni et al., 2002), estimates of vegetation biomass (Carvalhais et al., 2014), and soil organic carbon content (Hugelius et al., 2013a, b) (Table S3). We aggregated ∼130 000 1 km resolution MODIS LAI data monthly within each 1^∘ × 1^∘ pixel. Biomass based on remotely sensed forest biomass (Thurner et al., 2014) and upscaled GPP (Jung et al., 2011) covering the pan-Arctic domain was aggregated to 1^∘ resolution (Carvalhais et al., 2014). We used the NCSCD spatially explicit product (Hugelius et al., 2013a, b), which was generated from 1778 soil sample locations interpolated to a 1^∘ grid.

We apply the set-up described above to 3304 1^∘ × 1^∘ pixels (1686 in tundra; 1618 in taiga) using a monthly time step. Each pixel is treated independently without assuming a prior land cover or plant functional type, and we assume no spatial correlation between uncertainties in all pixels. In each 1^∘ × 1^∘ pixel, we applied the MHMCMC algorithm to determine the probability distribution of the optimal parameter set and initial conditions (x_i; Table S2) given observational constraints (O_i; LAI, SOC, and biomass; Table S3) using the same Bayesian inference approach described in Bloom et al. (2016):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }p\left(x_{\mathrm{i}}|O_{\mathrm{i}}\right)\propto p\left(x_{\mathrm{i}}\right)\times p\left(O_{\mathrm{i}}|x_{\mathrm{i}}\right).\end{array}$$

First, in the expression 1, p(x_i) represents the prior probability distribution of each DALEC2 parameter (x_i) and is expressed as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }p\left(x_{\mathrm{i}}\right)=p_{\mathrm{EDC}}\left(x_{\mathrm{i}}\right)\times e^{-\mathrm{0.5}{\left(\frac{\log \left(f_{\mathrm{auto}}\right)-\log \left(\mathrm{0.5}\right)}{\log \left(\mathrm{1.2}\right)}\right)}^{\mathrm{2}}}\times e^{-\mathrm{0.5}{\left(\frac{\log \left(C_{\mathrm{eff}}\right)-\log \left(\mathrm{17.5}\right)}{\log \left(\mathrm{1.2}\right)}\right)}^{\mathrm{2}}},\end{array}$$

where p_EDC(x_i) is the prior parameter probability according to the EDCs included in Table S2 and described in Bloom and Williams (2015). In addition, prior values for two parameters and their uncertainties (canopy efficiency, C_eff, and fraction of GPP respired, f_auto) are imposed with a log-normal distribution following Bloom et al. (2016) to be consistent with the global GPP range estimated in Beer et al. (2010) and f_auto ranges specified by DeLucia et al. (2007), respectively.

Second, p(O|x_i) from expression 1 represents the likelihood of x_i with respect to O_i, and it is calculated based on the ability of DALEC2 to reproduce (1) biomass (Carvalhais et al., 2014), (2) SOC (Hugelius et al., 2013a, b), and (3) MODIS LAI (Myneni et al., 2002). The reported uncertainty in biomass data from Thurner et al. (2014) was ±37 % at pixel scale. Because of undetermined errors related to tree cover thresholds used in the upscaling and to reflect unknown model structural error, we slightly inflate the error estimate and use a log-transform (1.5) of $\times /÷\mathrm{1.5}$ (i.e. $\times /÷\mathrm{1.5}$ spans 67 % of the expected error). We use the same proportional error for SOC. For MODIS LAI we inflate the proportional error further to log(2) based on well-reported biases in this product for evergreen forests (De Kauwe et al., 2011) and the estimated measurement and aggregation uncertainty for boreal forest LAI of 1 m² m⁻² reported by Goulden et al. (2011). The uncertainty assumptions in expression 3 are chosen due to a lack of better knowledge about the combined uncertainties arising from model representation errors and observation errors:

$$\begin{array}{ll}{\displaystyle }p& {\displaystyle }\left(O_{\mathrm{i}}|x_{\mathrm{i}}\right)=e^{-\mathrm{0.5}{\left(\frac{\log \left(O_{\mathrm{biomass}}\right)-\log \left(M_{\mathrm{biomass},\mathrm{0}}\right)}{\log \left(\mathrm{1.5}\right)}\right)}^{\mathrm{2}}}\\ \text{(3)}& {\displaystyle }& {\displaystyle }\times e^{-\mathrm{0.5}{\left(\frac{\log \left(O_{\mathrm{SOC}}\right)-\log \left(M_{\mathrm{SOC},\mathrm{0}}\right)}{\log \left(\mathrm{1.5}\right)}\right)}^{\mathrm{2}}}\times e^{-\mathrm{0.5}{\left(\frac{\log \left(O_{\mathrm{LAI},t}\right)-\log \left(M_{\mathrm{LAI},t}\right)}{\log \left(\mathrm{2}\right)}\right)}^{\mathrm{2}}}.\end{array}$$

For each 1^∘ × 1^∘ pixel we run three MHMCMC chains with 10⁷ accepted simulations each until convergence of at least two chains. We use 500 parameter sets sampled from the second half of each chain to describe the posterior distribution of parameter sets. We produce confidence intervals of terrestrial C fluxes and stocks from the selected parameter sets. In the following we report the highest confidence results (median; P₅₀) and the uncertainty represented by the 90 % confidence interval (5th percentile to 95th percentile – P₀₅, P₉₅). We calculate the transit time for C pools using the approach for non-steady-state pools described in Bloom et al. (2016), Sect. S3.

2.3 Model evaluation against independent in situ and pan-Arctic datasets

At the pan-Arctic scale, we compared CARDAMOM GPP with the FLUXCOM dataset from Jung et al. (2017). We also compared our CARDAMOM R_h with the global spatio-temporal distribution of soil respiration from Hashimoto et al. (2015) calculated by a climate-driven empirical model. To assess the degree of statistical agreement we calculated linear goodness of fit (slope, intercept, R²) between CARDAMOM and the two independent datasets and determined RMSE and bias from direct comparison on model–data residuals. The mapping includes stipples representing locations where the independent datasets are within CARDAMOM's 90 % confidence interval.

At a local scale, we compare CARDAMOM NEE and its partitioned components GPP and R_eco estimates with monthly aggregated values from the FLUXNET2015 sites. We selected eight sites (Belelli Marchesini et al., 2007; Bond-Lamberty et al., 2004; Goulden et al., 1996; Ikawa et al., 2015; Kutzbach et al., 2007; López-Blanco et al., 2017; Lund et al., 2012; Sari et al., 2017) located across sub- and high-Arctic latitudes, covering locations with different climatic conditions and dominating ecotypes (Table S4). For this evaluation, we compared the same years for both observations and CARDAMOM, and we selected data using the daytime method (Lasslop et al., 2010) due to the absence of a true night-time period during Arctic summers in some locations. Additionally, we selected a variable u_* threshold to identify insufficient turbulence wind conditions from year to year similar to López-Blanco et al. (2017). In this data–model comparison we included the median (P₅₀) ± the 90 % confidence interval (5th to 95th percentiles, (P₀₅, P₉₅)) including both random and u_* filtering uncertainty following the method described in Papale et al. (2006). Some of the sites lack wintertime measurements, and we filtered out data for months with less than 10 % observations. We performed a point-to-grid-cell comparison to assess the degree of agreement between each flux magnitude and seasonality calculating the statistics of linear fit (slope, intercept, R²) per flux and site between CARDAMOM and FLUXNET2015 datasets and determined RMSE and bias from model–data residuals comparison.

2.4 Benchmark of global vegetation models from ISI-MIP2a

We compared CARDAMOM analyses of pan-Arctic NPP, vegetation biomass carbon stocks (C_veg) and vegetation transit times (TT_veg) with six participating GVMs in the ISI-MIP2a comparison project (Akihiko et al., 2017). In this study we have considered DLEM (Tian et al., 2015), LPJmL (Schaphoff et al., 2013; Sitch et al., 2003), LPJ-GUESS (Smith et al., 2014), ORCHIDEE (Guimberteau et al., 2018), VEGAS (Zeng et al., 2005), and VISIT (Ito and Inatomi, 2012). The specific properties and degree of complexity of each ISI-MIP2a model are summarised in Table S5. The comparisons have been performed under the same spatial resolution as the CARDAMOM spatial resolution (1^∘ × 1^∘) for the 2000–2010 period. Also, the chosen GVMs from the ISI-MIP2a phase have their forcing based on ERA-Interim climate data, similar to the forcing used in CARDAMOM. We estimated the degree of agreement using the statistics of linear fit (slope, intercept, R², RMSE, and bias) per variable and model between CARDAMOM and GVMs but also their spatial variability including stipples where the GVM datasets are within CARDAMOM's 90 % confidence interval.

To understand the sources of errors in TT_veg calculations, we used CARDAMOM to calculate two hypothetical TT_veg (i.e. EXPERIMENT A TT_veg = ISI-MIP2a C_veg∕ CARDAMOM NPP and EXPERIMENT B TT_veg = CARDAMOM C_veg∕ ISI-MIP2a NPP) and then assessed the largest difference with CARDAMOM's CONTROL TT_veg. We estimated the hypothetical TT_veg for each pixel in each model and derived a pixel-wise measure of the contribution of biases in NPP and C_veg to biases in TT_veg by overlapping their distribution functions.

3 Results

3.1 Pan-Arctic retrievals of C cycle

Overall, we found that the pan-Arctic region (Fig. 1 and Table 1) acted as a small sink of C (area-weighted P₅₀) over the 2000–2015 period with an average of $-\mathrm{67.4}(-\mathrm{286.7}$, 1159.9) g C m⁻² yr⁻¹, P₅₀ (P₀₅, P₉₅), although the 90 % confidence intervals remain large (and so the region could be a source of C). Tundra region NEE was estimated at $-\mathrm{14.9}(-\mathrm{163.4}$, 1116.1) g C m⁻² yr⁻¹, a weaker sink compared to taiga regions, which were established at $-\mathrm{110.4}(-\mathrm{387.7}$, 1195.8) g C m⁻² yr⁻¹. The photosynthetic inputs exceeded the respiratory outputs (GPP  > R_eco; Table 1), although the much larger uncertainties stemming from R_eco, and more specifically from R_h, compared with GPP complicate the net C sink–source estimate beyond the median's average ensembles. In the pan-Arctic region approximately half of GPP is autotrophically respired, resulting in an NPP of 290.3 (196.4, 410.7) g C m⁻² yr⁻¹. Carbon use efficiency (NPP∕GPP) averages 0.51 (0.46, 0.55), and varied marginally across tundra 0.51 (0.46, 0.54) and taiga 0.52 (0.46, 0.56). Despite these apparent small variations, tundra photosynthesised and respired (respectively, 327.2 (236.8, 463.3) and 310.0 (124.3, 1536.8) g C m⁻² yr⁻¹) approximately half as much as the Taiga region (759.8 (584.1, 967.9) and 635.3 (285.3, 2114.0) g C m⁻² yr⁻¹).

Figure 1 Schematic diagram of the terrestrial C processes modelled in CARDAMOM for the pan-Arctic (black values), tundra (yellow values), and taiga (green values) domains. The values characterise the median for the 2000–2015 period and the parentheses delimit the 90 % confidence interval. C processes represented include flows for C fluxes in white (NEE, net ecosystem exchange; GPP, gross primary production; NPP, net primary production; R_eco, ecosystem respiration; R_a, autotrophic respiration; R_h, heterotrophic respiration), C allocation in dark green (to labile, leaf, stem and root), and C turnover in cyan (from leaf, wood, roots, and litter). C stocks are represented in dark blue boxes (labile, leaf, stem, root, litter, and SOM, soil organic matter) and aggregated into photosynthetic (C_photo = leaf + labile), vegetation (C_veg = leaf + labile + wood + roots), soil (C_dom = litter + SOM), and total (C_tot = C_photo + C_veg + C_dom) C stocks in red boxes. Analogy, transit times (TT) are also aggregated into photosynthetic (TT_photo = leaf + labile), vegetation (TT_veg = leaf + labile + wood + roots), soil (TT_dom = litter + SOM), and total (TT_tot = TT_photo + TT_veg + TT_dom) C transit times.

Table 1 Multi-year (2000–2015) annual average of main ecosystem C fluxes (NEE, GPP, NPP, R_eco, R_a, R_h; g C m⁻² yr⁻¹), C stocks (C_photo, C_veg, C_dom, C_tot; kg C m⁻²) and transit times (TT_photo, TT_veg, TT_dom, TT_tot; years) for the pan-Arctic, tundra (non-forested), and taiga (forested) domain. The averages contain the median in bold (50th percentile) and the uncertainty estimations across the 90 % confidence range between the 5th and 95th percentiles assuming no spatial correlation between uncertainties in all pixels. We assume spatial correlation between pixels: P₅₀, P₀₅, and P₉₅ represent the area-weighted aggregate of all pixels' median, P₀₅, and P₉₅.

The total size of the pan-Arctic soil C stock (C_dom) averaged 24.4 (10.3, 47.5) kg C m⁻², 16-fold greater than the vegetation C stock (C_veg), 1.5 (0.5, 5.8) kg C m⁻². The soil C stock (fresh litter and SOM) is dominated by C_som, accounting for the 99 % which also dominates the total terrestrial C stock in the pan-Arctic. Among the living C stocks, 93 % of C (88 % in tundra and 90 % in taiga) is allocated to the structural stocks (wood and roots; 1.4 (0.4, 5.6) kg C m⁻²) compared to 7 % (12 % in tundra and 10 % in taiga) to the photosynthetic stock (leaves and labile; 0.1 (0.1, 0.2) kg C m⁻²). On average, the total ecosystem C stock is 26.3 (11.8, 51.0) kg C m⁻² in the pan-Arctic region, with slightly lower stocks in tundra (24.6 (10.8, 50.6) kg C m⁻²) than taiga (27.7 (12.7, 51.2) kg C m⁻²). In general, the taiga region holds on average ∼100 % more photosynthetic tissues, ∼160 % more structural tissue, and ∼9 % more soil C stocks than tundra. In other words, taiga holds ∼12 % more total C than tundra. The greater living stock of C in taiga (2.1 (0.8, 5.1) kg C m⁻²) than tundra (0.8 (0.3, 6.8) kg C m⁻²) means that the relative size of R_a and R_h in the two regions differs. Thus in tundra R_a accounts for 51 % of total ecosystem respiration, while in taiga this fraction is 57 %. R_a is 4 % larger than R_h in tundra but 24 % greater in taiga, reflecting the greater rates of C cycling in taiga. Uncertainties in estimates of soil C stock are notably higher than for living C stocks, highlighting the lack of observational and mechanistic constraint on heterotrophic respiration.

The global mean C transit time is 1.3 (0.8, 2.1) years in leaves and labile plant tissue (TT_photo), 4.5 (1.7, 15.7) years in stems and roots (TT_veg), and 120.5 (9.8, 822.6) years in litter and SOM (TT_dom). The total C transit time (TT_tot) (133.1 (11.5, 1013.6) years) is clearly dominated by the soil C stock, highlighting the very long periods of times that C persists in Arctic soils. CARDAMOM calculated 62 % longer TT_dom in tundra compared to taiga, likely linked to lower temperatures, but uncertainties are large due to the limitations of data constraints.

Figure 2 Original soil organic carbon (SOC; Hugelius et al., 2013a), biomass (Carvalhais et al., 2014), and leaf area index (LAI; Myneni et al., 2002) datasets used in the data-assimilation process within the CARDAMOM framework (a–c), assimilated SOC, biomass, and LAI integrated into CARDAMOM (d–f) and their respective goodness-of-fit statistics between original and assimilated datasets (g–i). The error bars represent the 90 % confidence interval of the assimilated variable in CARDAMOM.

3.2 Data assimilation and uncertainty reduction

The CARDAMOM framework generated an analysis broadly consistent with the combination of SOC, biomass, and LAI in each grid cell (Fig. 2) and the errors assigned to these data products. The agreement for the SOC dataset by Hugelius et al. (2013a) is a 1:1 relationship (R²=1.0; RMSE = 0.95 kg C m⁻²), reflecting a straightforward model parameterisation. The biomass product from Carvalhais et al. (2014) was well correlated (R²=0.97; RMSE = 0.46 kg C m⁻²), but CARDAMOM was consistently biased ∼28 % low. MODIS LAI data were also well correlated ($R^{\mathrm{2}}=\mathrm{0.79};$ RMSE = 0.42 kg C m⁻²) but ∼28 % higher than CARDAMOM analyses. These biases (Fig. 2) likely arise due to a low estimate in the photosynthesis model (ACM) used in CARDAMOM which propagates through the C cycle. CARDAMOM balances uncertainty in data products and the models (ACM photosynthesis model and DALEC2) to generate a weighted analysis typical of Bayesian approaches. The CARDAMOM analysis 90 % confidence interval (CI) includes the 1:1 line for biomass and LAI (Fig. 2), indicating that the likelihoods on C cycle analyses include the expected value of the observations.

Figure 3 Original gross primary productivity (GPP; Jung et al., 2017) and heterotrophic respiration (R_h; Hashimoto et al., 2015) datasets used in the data validation process (a, b), estimated GPP and R_h by CARDAMOM (c, d), and their respective goodness-of-fit statistics between original and assimilated datasets (e, f). Stippling indicates locations where the independent datasets are within CARDAMOM's 5th and 95th percentiles.

The degree to which posterior distributions were constrained from the prior distributions in each of the 17 model parameters and 6 initial stock sizes (Table S2) varied considerably depending on the parameters in question and their related processes (Table 2 and Fig. S2). The 90 % CI posterior range of foliar, wood, labile, and SOM C stocks (C_foliar, C_wood, C_labile, and C_som) as well as parameters such as allocation to foliage (f_fol) and lifespan (L) were considerably reduced (>80 % uncertainty reduction compared to priors) most likely controlled by the information on LAI, biomass, and SOC constraints. Contrarily, parameters that have not been regulated in any way in the MHMCMC algorithm, i.e. turnover processes such as litter mineralisation (MR_litter), root turnover (TOR_roots), wood turnover (TOR_wood), decomposition rates (D_rate), and initial C stock such as litter (C_litter) were found to be poorly constrained (<20 % uncertainty reduction). Overall, the uncertainty reduction classified by processes and ranked from most to least constrained estimated a 71 % reduction for C stocks, a 67 % reduction for C allocation, 59 % for plant phenology, and 31 % for C turnover related parameters. Although there are no substantial differences between tundra and taiga, C_roots was better constrained in tundra regions (42 %), while leaf onset day (B_day), leaf fall day (F_day), and leaf fall duration (L_f) were better constrained in taiga regions (>18 % or more).

Table 2 Parameter uncertainty reduction in percentage ranked from least (red) to most (blue) constrained in the pan-Arctic, tundra, and taiga domains. The reduction percentage is calculated based on the difference between the 90 % CI prior range and the 90 % CI posterior range.

3.3 Independent evaluation: from global to local scale

We compared our estimates of GPP and R_h with independent datasets to evaluate the model performance (Fig. 3). We found GPP to be well correlated (R²=0.81; RMSE = 0.43 kg C m⁻²) but biased lower (∼53 %) compared to Jung et al.'s (2017) GPP estimates. There are in general very few pixels where the FLUXCOM product falls within CARDAMOM's 90 % confidence interval. Additionally, the R_h product from Hashimoto et al. (2015) is less consistent with our estimates (R²=0.40; RMSE = 0.09 kg C m⁻²), presenting a tendency towards lower values in tundra pixels and higher values in taiga pixels. The spatial variability of R_h is considerably smaller in Hashimoto et al. (2015) compared to our CARDAMOM estimates. R_h falls within the 90 % confidence interval of CARDAMOM in most of the pan-Arctic region due to the fact that the R_h uncertainties are significant (Fig. 3). This finding confirms the uncertainties previously noted in modelled respiratory processes (Table 1) where the upper P₉₅ in R_h dominated NEE's uncertainties but also the soil C stocks and transit times.

Figure 4 Monthly aggregated seasonal variability of observed (FLUXNET2015) and modelled (CARDAMOM) C fluxes (NEE, net ecosystem exchange; GPP, gross primary production; R_eco, ecosystem respiration) across eight low- and high-Arctic sites (Hakasia, Kobbefjord, Manitoba, Poker Flat, Samoylov, Tiksi, UCI-1998, and Zackenberg). Each of these sites, located in different countries (RU – Russia; GL – Greenland; CA – Canada; US – United States) feature different meteorological conditions and vegetation types (Table S4). Uncertainties represent the 25th and 50th percentiles (darker shade) and the 5th and 95th percentiles (lighter shade) of both field observations and the CARDAMOM framework.

For comparison with direct ground observations from the FLUXNET2015 dataset, we report here monthly aggregated $P_{\mathrm{50}}\pm P_{\mathrm{05}-\mathrm{95}}$ estimates of NEE, GPP, and R_eco to show timing and magnitudes but also to diagnose whether CARDAMOM is in general agreement with flux tower data. Overall, CARDAMOM performed well in simulating observed NEE ($R^{\mathrm{2}}=\mathrm{0.66};$ RMSE = 0.51 g C m⁻² per month; bias = 0.16 g C m⁻² per month), GPP (R²=0.85; RMSE = 0.89 g C m⁻² per month; bias = 0.5 g C m⁻² per month) and R_eco (R²=0.82; RMSE = 0.63 g C m⁻² per month; bias = 0.35 g C m⁻² per month) across eight sub-Arctic and high-Arctic sites from the FLUXNET2015 dataset (Fig. 4; Table S6). CARDAMOM NEE is ∼25 % lower than FLUXNET2015, while GPP and R_eco are ∼30 % and ∼10 % higher, respectively. This mismatch is important in the context of the FLUXCOM GPP upscaling, 50 % higher than CARDAMOM GPP. At some sites such as Hakasia, Samoylov, Poker Flat, and Manitoba (NEE R²=0.73; GPP R²=0.92; R_eco R²=0.88), CARDAMOM better matches the seasonality and the magnitude of the C fluxes than the rest, i.e. Tiksi, Kobbefjord, Zackenberg, and UCI-1998 (NEE R²=0.58; GPP R²=0.67; R_eco R²=0.67). In general, CARDAMOM captured the beginning and the end of the growing season well (Fig. 4), although the assimilation system has some bias due to (1) a difference in timing (e.g. earlier shifts of peak of the growing season in Manitoba GPP and R_eco and earlier end of the growing season in Poker Flat NEE) and (2) differences in flux magnitudes (such as in Hakasia GPP and R_eco and Kobbefjord NEE).

3.4 Benchmarking ISI-MIP2a models with CARDAMOM

We used our highest confidence retrievals of NPP, C_veg, and TT_veg (i.e. retrievals including assimilated LAI, biomass, and SOC) to benchmark the performance of the GVMs from the ISI-MIP2a project. In this assessment we compared not only their spatial variability across the pan-Arctic, tundra, and taiga region (Fig. 5) but also the degree of agreement between their mean model ensemble within the 90 % confidence interval of our assimilation framework (Fig. 6, Table 3). NPP estimates (RMSE = 0.1 kg C m⁻² yr⁻¹; R²=0.44) are in better agreement than C_veg (RMSE = 1.8 kg C m⁻²; R²=0.22) and TT_veg (RMSE = 4.1 years; R²=0.12). The assessed GVMs estimated on average 8 % lower NPP, 16 % higher C_veg, and 22 % longer TT_veg than CARDAMOM across the entire pan-Arctic domain (Figs. 5 and 6) on average. Thus, at regional aggregation CARDAMOM analyses agreed more closely with ISI-MIP2a models than with FLUXCOM (51 % difference) and with the Carvalhais et al. (2014) biomass data (28 % bias).

Figure 5 Central tendency and variability of NPP (net primary production), C_veg (vegetation C stock), TT_veg (vegetation transit time) estimated by CARDAMOM (orange), and ISI-MIP2a models (grey) in the pan-Arctic, tundra, and taiga regions. The box–whisker plots comprise the estimations between the 5th and 95th percentiles, and the box encompasses the 25th to 75th percentiles. The line in each box marks the median of studied variables in each region.

Figure 6 NPP (net primary production), C_veg (vegetation C stock), and TT_veg (vegetation transit time) ratios between ISI-MIP2a model ensembles (DLEM, LPJmL, LPJ-GUESS, ORCHIDEE, VEGAS, and VISIT) and CARDAMOM. Stippling indicates locations where the ISI-MIP2a model mean is within CARDAMOM's 5th and 95th percentiles.

Table 3 Statistics of linear fit between the CARDAMOM framework (independent) and the ISI-MIP2a models (dependent) per individual model and per NPP (net primary production; kg C m⁻² yr⁻¹), C_veg (vegetation C stock; kg C m⁻²), and TT_veg (vegetation transit time; years). The units for RMSE and bias are kg C m⁻² yr⁻¹ in NPP, kg C m⁻² yr⁻¹ in C_veg, and years in TT_veg.

The poor spatial agreement regarding TT_veg between CARDAMOM and ISI-MIP2a (Table 3) is indicative of uncertainties in the internal C dynamics of these models. For instance, the slopes in Table 3 are steep and the R² are poor – so there is substantial disagreement in the spatial pattern, not just a large bias. For ISI-MIP2a comparison R² values ranged from 0.03 to 0.52 for NPP, from 0.00 to 0.31 for C_veg, and from 0.00 to 0.24 for TT_veg. Spatially, the stippling in Fig. 6 indicates areas where the GVMs are within the 90 % CI of CARDAMOM; agreement is best over the taiga domain rather than in tundra for TT_veg. The benchmark area of consistency (stippling) is more extensive for C_veg and TT_veg than for NPP. Thus, while there is a stronger spatial correlation for NPP between CARDAMOM and GVMs (Table 3), this is a clearer bias for NPP. Some models (LPJ-GUESS and ORCHIDEE) systematically calculate lower values in all the assessed variables, while others (LPJmL and VISIT) calculate higher estimates. The models in closer agreement with CARDAMOM were DLEM (5 % difference) and LPJ-GUESS (17 %), while VEGAS (44 %) and ORCHIDEE (56 %) were the models with larger discrepancies (Table 3; Figs. 5 and 6).

The attribution analysis to identify the origin of bias from ISI-MIP2a models indicated a joint split between NPP and C_veg for TT_veg error simulated in GVMs (Fig. 7). The distribution of the differences relative to CARDAMOM revealed that the higher error (i.e. the lower overlapped area, and by extension the largest contributor to TT_veg biases) comes from ISI-MIP2a NPP with a 69 % agreement in the distribution, while C_veg agrees to 72 %. In fact, the TT_veg R² for each model (Table 3) is very close to the product of the NPP R² and C_veg R² for that model, i.e. the uncertainty in the TT_veg is a direct interaction of NPP and C_veg uncertainty (R² of the correlation: 0.71). This finding supports Fig. 6, which shows TT_veg error derives equally from both NPP and C_veg.

Figure 7 Distribution functions derived from the attribution analysis used to estimate the origin of vegetation transit time (TT_veg) bias from ISI-MIP2a models. The CONTROL TT (grey) includes both biomass (C_veg) and net primary production (NPP) estimated by CARDAMOM. EXPERIMENT A TT (dark red) incorporates C_veg from ISI-MIP2a and NPP from CARDAMOM, while EXPERIMENT B TT (dark green) includes NPP from ISI-MIP2a and C_veg from CARDAMOM. The lower the overlapped area between control and experimental TT, the larger the contribution for TT biases. For readability purposes, the scale on the x axis is limited to 20 years.

4 Discussion

4.1 Pan-Arctic retrievals of C cycle

The CARDAMOM framework has been used to evaluate the terrestrial pan-Arctic C cycle in tundra and taiga at coarse spatio-temporal scale (at monthly and annual time steps for the 2000–2015 period and at 1^∘ × 1^∘ grid cells). Overall, we found that the pan-Arctic region was most likely a consistent sink of C (weaker in tundra and stronger in taiga), although the large uncertainties derived from respiratory processes (Table 1) strongly increase the 90 % confidence interval uncertainty. We estimate that tundra experienced 62 % longer transit times in litter and SOM C stocks than taiga ecosystems. Further, the contribution of R_a and R_h to total ecosystem respiration was similar in tundra (51 %, 49 %) but dominated by R_a in taiga (57 % compared to 43 %).

CARDAMOM retrievals are consistent with outcomes from relevant papers such as the (i) C flux observations and model estimates reported in McGuire et al. (2012), (ii) C stocks and transit times described by Carvalhais et al. (2014), and (iii) NPP, C stocks, and turnover rates stated in Thurner et al. (2017).

• i.

  The CARDAMOM NEE estimates reported in this study for the tundra domain are inside the variability comparison of values compiled by McGuire et al. (2012) considering field observation, regional-process-based models, global-process-based models, and inversion models. The authors reported that Arctic tundra was a sink of CO₂ of −150 Tg C yr⁻¹ (SD = 45.9) across the 2000–2006 period over an area of 9.16×10⁶  m². Here, CARDAMOM NEE estimated −129 Tg C yr⁻¹ over an area of 8.1×10⁶ km² for the same period. This exhaustive assessment of the C balance in Arctic tundra included approximately 250 estimates using the chamber and eddy covariance method from 120 published papers (McGuire et al., 2012; Supplement 1) with an area-weighted mean of means of −202 Tg C yr⁻¹. The regional models, including runs from LPJ-WHyMe (Wania et al., 2009a, b), ORCHIDEE (Koven et al., 2011), TEM6 (McGuire et al., 2010), and the TCF model (Kimball et al., 2009), reported an NEE of −187 Tg C yr⁻¹ and GPP, NPP, R_a, and R_h of 350, 199, 151, and 182 g C m⁻² yr⁻¹, respectively. GVMs applications such as CLM4C (Lawrence et al., 2011), CLM4CN (Thornton et al., 2009), Hyland (Levy et al., 2004), LPJ (Sitch et al., 2003), LPJ-GUESS (Smith et al., 2001), O-CN (Zaehle and Friend, 2010), SDGVM (Woodward et al., 1995), and TRIFFID (Cox, 2001) estimated an NEE of −93 Tg C yr⁻¹ and GPP, NPP, R_a, and R_h of 272, 162, 83, and 144 g C m⁻² yr⁻¹. For the same period, CARDAMOM has estimated 330, 167, 160, and 154 g C m⁻² yr⁻¹, respectively, for the same gross C fluxes.

• ii.

  Carvalhais et al. (2014) estimated a total ecosystem carbon (C_tot) of 20.5 (8.0, 52.5) kg C m⁻² for tundra and 24.8 (15.2, 58.0) kg C m⁻² for taiga, while values from CARDAMOM were 24.6 (10.8, 50.6) kg C m⁻² for tundra and 27.7 (12.7, 51.2) kg C m⁻² in taiga (Fig. 5; Table 1) for the same area. Thus, Carvalhais et al.'s (2014) C_tot product stored 20 % and 12 % less carbon in tundra and taiga, respectively, than CARDAMOM. Overall, CARDAMOM calculated 20 % and 6 % longer transit times for tundra and taiga, respectively, with average values of 80.8 (21.8, 195.2) years in tundra and 51.2 (22.1, 109.3) years in taiga (Table 1) compared to the 64.4 (25.7, 259.8) years in tundra and 48.2 (111.6, 24.9) years in taiga in Carvalhais et al. (2014). These numbers have been retrieved from the same biome classification and they include the 90 % confidence interval of the assessed spatial variability. Also, we applied a correction factor of TT_gpp = TT_npp × (1 − fraction of GPP respired) to be comparable with Carvalhais et al. (2014) TT. Both datasets agree on the fact that at high (cold) latitudes, first tundra and secondly taiga have the longest transit times on the entire globe (Bloom et al., 2016; Carvalhais et al., 2014).

• iii.

  A recent study from Thurner et al. (2017) assessed temperate and taiga-related TTs, presenting a 5-year average NPP dataset applying both MODIS (Running et al., 2004; Zhao et al., 2005) and BETHY/DLR (Tum et al., 2016) products and an innovative biomass product (Thurner et al., 2014) accounting for both forest and non-forest vegetation. Our estimate of TT_veg for the exact same period is 5.3 (1.9, 18.2) years, compared to Thurner et al.'s (2017) TT, 8.2 (5.5, 11.5) years using MODIS, and 6.5 (4.2, 8.7) years using BETHY/DLR. A note of caution here: the number reported by the authors are turnover rates, which are converted to transit times by just applying the inverse of turnover rates (TT${}_{\mathrm{veg}}=\mathrm{1}/$turnover rates). Additionally, their NPP estimates, 0.35 and 0.45 kg C m⁻² yr⁻¹ from both MODIS and BETHY/DLR, are only 5 % more productive on average than the CARDAMOM NPP estimate (0.4 (0.3, 0.6) kg C m⁻² yr⁻¹). The biomass derived from Thurner et al. (2014) (3.0±1.1 kg C m⁻²) is ∼30 % lower than CARDAMOM C_veg (2.2 (1.1, 5.0) kg C m⁻²), calculated for the same period and for the same taiga domain.

In general, we found a reasonable agreement between CARDAMOM and assimilated and independent data at pan-Arctic scale. CARDAMOM retrievals of assimilated data are in good agreement with the SOC (Fig. 2). The simulation of TT_dom is weakly constrained (Table 1) – our analysis adjusts TT to match mapped stocks, hence the strong match of modelled to mapped SOC. So, independent data on TT_dom data (e.g. ¹⁴C) are required across the pan-Arctic region to provide a stronger constraint on process parameters and reduce the very broad confidence intervals of CARDAMOM analyses. The low bias in mean estimates of LAI and biomass (Fig. 2) likely relates to the strong prior on photosynthesis estimates from the ACM model, which lacks a temperature acclimation for high latitudes in this implementation. However, the uncertainty in the biomass and LAI analyses spans the magnitude of the bias. So, CARDAMOM generates some parameters sets that are consistent with observations. CARDAMOM produces analyses that reproduce the pattern of LAI, GPP, biomass, and SOC (Figs. 2 and 3) – this demonstrates that the DALEC model structure can be calibrated to simulate the links between these variables as a function of mass balance constraints and realistic process interactions and climate sensitivities.

There are clear biases in CARDAMOM analyses compared to independent global upscaled GPP (Jung et al., 2017) and R_h products (Hashimoto et al., 2015) (Fig. 3). However, CARDAMOM resolves the spatial pattern in GPP effectively, while the spatial mismatch in R_h estimates is marked (Fig. 3), echoing the large uncertainty found in NEE (Fig. 1, Table 1). One difference with Hashimoto et al.'s (2015) R_h model is the lack of moisture limitation on respiration in CARDAMOM. Conversely, GPP is relatively well-constrained in space through the assimilation of LAI and a prior for productivity (Bloom et al., 2016), although an important mismatch has been found: CARDAMOM GPP is 50 % lower than FLUXCOM but 30 % higher than FLUXNET2015 EC data.

The agreement between CARDAMOM analyses and EC data is high given the scale difference. A direct point-to-grid cell comparison with local observations derived from the FLUXNET2015 dataset (Fig. 4, Table S6) is challenging and always difficult. CARDAMOM outputs covers 1^∘ × 1^∘ grid cells, whereas local eddy covariance flux measurements are of the order of 1–10 ha. Thus, for observational sites located in areas with complex terrain, such as Kobbefjord in coastal Greenland, the agreement can be expected to be low. For inland forest sites, such as Poker Flat in Alaska, there may be less differences in vegetation characteristics and local climatology between the local-scale measurement footprint and the corresponding CARDAMOM grid cell. This scaling issue is likely to have a larger impact on flux magnitudes compared with seasonal dynamics. In general, CARDAMOM captured the seasonal dynamics in NEE, GPP, and R_eco well (Fig. 4, Table S6), although the monthly model time step does reduce skill in shoulder seasons. There was a consistent timing mismatch in early season flux increase, where CARDAMOM predicts earlier growing season onset compared with observations. This is likely due to the impact of snow cover, which is not explicitly included in the CARDAMOM framework.

For a further independent evaluation of CARDAMOM outputs, we compare the tundra and boreal estimates to plot scale flux and stock information. For tundra, Street et al. (2012) calculate growing season GPP estimates of 263–380 g C m⁻² for Empetrum nigrum communities and 295–386 g C m⁻² for Betula nana communities, which is consistent with the ranges in Fig. 1 for tundra. Biomass stocks for Arctic tundra recorded in the Arctic LTER (Long Term Ecological Research) at Toolik Lake range from 105 to 1160 g C m⁻² (Hobbie and Kling, 2014), which are consistent with the estimates from CARDAMOM, albeit at the lower end of the model estimates. For boreal forests, Goulden et al. (2011) report annual GPP estimates across a chronosequence of stands, and thus a variation across canopy densities, which varied from 450 to 720 g C m⁻² yr⁻¹. These data are consistent with the span of GPP in CARDAMOM (Fig. 1), again best matching the lower end of the model estimates. For the same study, the vegetation C stock estimates varied from 100 to 5000 g C m⁻², consistent with CARDAMOM and with measurements of 10 to 40-year old boreal stands best matching the CARDAMOM median estimate of ∼1500 g C m⁻². We conclude from comparisons with site data that CARDAMOM analyses are broadly consistent, with some tendency for CARDAMOM to have a high bias. This comparison is similar to the FLUXNET2015 evaluation of CARDAMOM. But it conflicts with the estimation of low bias from the comparison of CARDAMOM with FLUXCOM GPP and Carvalhais et al. (2014) biomass stock maps. It is possible that the scale differences between site-level products and landscape estimates is confusing these comparisons, but there is clearly a need to understand these inconsistencies in C cycle estimates better.

4.2 CARDAMOM as a model benchmarking tool

An ideal benchmarking tool for GVMs would compare model state variables and fluxes with multiple, independent, unbiased, error-characterised measurements collected repeatedly at the same temporal and spatial resolution. Of course direct measurements of key C cycle variables like these are not available. Even at FLUXNET sites GPP and R_eco must be inferred, and NEE data are often gap-filled. Satellite data can provide continuous fields but do not directly measure ecological variables like biomass or LAI, so calibrated models are required to generate ecological products. Atmospheric conditions can introduce biases and data gaps into optical data that are poorly quantified. Upscaling of FLUXNET data requires other spatial data, e.g. MODIS LAI, which challenge the characterisation of error and generates complex hybrid products. We suggest that CARDAMOM provides some of the requirements of the ideal benchmark system – an error-characterised, complete analysis of the C cycle that is based on a range of observational products. CARDAMOM includes its own C cycle model; this has the advantage of evaluating the observational data for consistency (e.g. with mass balance), propagating error across the C cycle, and generating internal model variables such as TT. Further, the model is of intermediate complexity and independent of the benchmarked models.

Using CARDAMOM as a benchmarking tool for six GVMs we found disagreements that varied among models for spatial estimates of NPP, C_veg, and TT_veg across the pan-Arctic (Fig. 6) in comparison with CARDAMOM confidence intervals. GVM NPP estimates had a higher correlation than TT_veg and C_veg with CARDAMOM analyses (Table 3), but because CARDAMOM confidence intervals on NPP were relatively narrow (Fig. 1), the benchmarking scores from GVM NPP were relatively poor (Fig. 6). Consequently, we used CARDAMOM to calculate the relative contribution of productivity and biomass to the transit times bias by applying a simple attribution analysis (Fig. 7). We concluded that the largest bias to transit times originated not from a deficient understanding of one single component, but from an equal combination of both productivity and biomass errors together. Therefore, this study partially agrees with previous studies (Friend et al., 2014; Nishina et al., 2014; Thurner et al., 2017) highlighting the deficient representation of transit times or turnover dynamics, but we further suggest that global vegetation and Earth system modellers need to focus on the productivity and vegetation C stocks dynamics to improve inner C dynamics. A major challenge for GVMs is the spin-up problem (Exbrayat et al., 2014). GVMs need to find a way to ensure that the spin-up process produces biomass estimates consistent with the growing availability of biomass maps from Earth observations. CARDAMOM solves this problem by avoiding spin-up. Its fast run time allows the biomass maps to act as a constraint on the probability distribution of model parameters. There may be opportunities to use CARDAMOM-style approaches to assist the GVM community address this problem.

4.3 Outlook

Although CARDAMOM estimates for pan-Arctic C cycling are in moderately good agreement with observations and data constraints, we have not included important components controlling ecosystem processes that could potentially improve our understanding on C feedbacks, with an emphasis on high-latitude ecosystems. For example, thaw and the release of permafrost C is not represented in CARDAMOM, but the influence on vegetation dynamics, permafrost degradation, and soil respiration is critical at high latitudes (Koven et al., 2015; Parazoo et al., 2018). Also, Koven et al. (2017) have shown that soil thermal regimes are key to resolving the long-term vulnerability of soil C. Moreover, we have not characterised snow dynamics nor the insulating effect of snow affecting respiratory losses across wintertime periods (López-Blanco et al., 2018). Further, methane emissions, another important contributor to total C budget (Mastepanov et al., 2008; Zona et al., 2016), were neglected from this modelling exercise since it is not easy to model due to its three complex transport mechanisms (Walter et al., 2001).

However, our approach to use an intermediate-complexity model has the strong advantage of allowing very large (10⁷) model ensembles per pixel and thus a thorough exploration of model–parameter interactions that is not feasible with typical GVMs. Other viable options include using emulators (Fer et al., 2018) and particle filters (Arulampalam et al., 2002), but MCMC methods provide the most detailed description of error distributions. There remains a strong argument to utilise intermediate-complexity models like DALEC2 to evaluate the minimum level of detail required to represent ecosystem processes consistent with local observations and to allow testing of alternate model structures. Assimilating further data products, for instance patterns in soil hydrology and snow states across the pan-Arctic from Earth observation, could provide useful information on spatio-temporal controls on soil activity and microbial metabolism to constrain below-ground processes. This information would need to be tied to process-level information on SOM turnover generated from experimental studies and included in updated versions of DALEC. Thus, more field observations are crucial across the pan-Arctic, specifically on the decomposition and TT of SOC (He et al., 2016) and respiratory processes such as the partitioning of R_eco into R_a and R_h (Hobbie et al., 2000; McGuire et al., 2000) across the growing season and also during wintertime (Commane et al., 2017; Zona et al., 2016).

Our approach has used estimated observation error and inflated this to include unknown errors associated with model process representation. We currently lack any better knowledge of the combined uncertainties arising from model representation errors and observation errors. We acknowledge that all models are an imperfect representation of C dynamics, which generates irreconcilable model–data errors due to the inherent assumptions in model structure. Future analyses should investigate model structural error, using, for example, error-explicit Bayesian approaches (Xu et al., 2017), or comparing the likelihoods of alternate model structures of varying complexity. Using multiple sources of data, we have highlighted systematic errors in the model at a landscape scale (Figs. 2 and 3) for LAI, GPP, and biomass. However, these biases are not consistent for site-scale evaluations. Thus, a next step would be to include explicitly both random and systematic process errors for C fluxes in the data assimilation. These errors could be determined from field-scale evaluation of model process representation (Table 2) using, e.g., FLUXNET2015 data. We also need to understand better the error associated with the landscape heterogeneity of C stocks and fluxes to upscale from flux tower observations or direct measurements of LAI to landscape pixel. This could be achieved by constructing robust observation error models (Dietze, 2017) from field to pixel scale for, e.g., GPP, LAI, and foliar N. The evaluation of the sensitivity of C cycling DA analyses to observation error has shown relatively low sensitivity to data gaps and random error on net ecosystem flux data (Hill et al., 2012), but further analyses of error sensitivity are required for multiple streams of stock data.

5 Conclusions

The Arctic is experiencing rapid environmental changes, which will influence the global C cycle. Using a data-assimilation framework we have evaluated the current state of key C flux, stocks, and transit times for the pan-Arctic region for 2000–2015. We found that the pan-Arctic was a likely sink of C, weaker in tundra and stronger in taiga, but uncertainties around the respiration losses are still large, and so the region could be a source of C. Comparisons with global- and local-scale datasets demonstrate the capabilities of CARDAMOM for analysing the C cycle in the Arctic domain. CARDAMOM is a data-constrained and data-integrated analysis, evaluated for internal consistency, and is therefore a good candidate to benchmark performance of global vegetation models. We conclude that a GVM bias found in transit time of vegetation C is the result of a joint combination of uncertainties from productivity processes and biomass in GVMs, and thus these are a major component of error in their forecasts. While spatial patterns in GVM predictions of NPP are reasonable, particularly in taiga, they have significant biases against the CARDAMOM benchmark. Improved mapping of vegetation and soil C stocks and change over time is required for better analytical constraint. Moreover, future work is required on assimilating data on soil hydrology, permafrost, and snow dynamics to improve accuracy and decrease uncertainties in belowground processes. This work establishes the baseline for further process-based ecological analyses using the CARDAMOM data-assimilation system as a technique to constrain the pan-Arctic C cycle.