Microbial decomposition processes and vulnerable arctic soil organic carbon in the 21st century

Various levels of representations of biogeochemical processes in current biogeochemistry models contribute to a large uncertainty in carbon budget quantification. Here, we present an uncertainty analysis with a process-based biogeochemistry model, the Terrestrial Ecosystem Model (TEM), into which detailed microbial mechanisms were incorporated. Ensemble regional simulations with the new model (MIC-TEM) estimated that the carbon budget of the arctic ecosystems is 76.0±114.8 Pg C during the 20th century, i.e., $-\mathrm{3.1}\pm \mathrm{61.7}$ Pg C under the RCP 2.6 scenario and 94.7±46 Pg C under the RCP 8.5 scenario during the 21st century. Positive values indicate the regional carbon sink while negative values are a source to the atmosphere. Compared to the estimates using a simpler soil decomposition algorithm in TEM, the new model estimated that the arctic terrestrial ecosystems stored 12 Pg less carbon over the 20th century, i.e., 19 and 30 Pg C less under the RCP 8.5 and RCP 2.6 scenarios, respectively, during the 21st century. When soil carbon within depths of 30, 100, and 300 cm was considered as initial carbon in the 21st century simulations, the region was estimated to accumulate 65.4, 88.6, and 109.8 Pg C, respectively, under the RCP 8.5 scenario. In contrast, under the RCP 2.6 scenario, the region lost 0.7, 2.2, and 3 Pg C, respectively, to the atmosphere. We conclude that the future regional carbon budget evaluation largely depends on whether or not adequate microbial activities are represented in earth system models and on the sizes of soil carbon considered in model simulations.

1 Introduction

Northern high-latitude soils and permafrost contain more than 1600 Pg carbon (Tarnocai et al., 2009). Climate over this region has warmed in recent decades (Serreze and Francis, 2006) and the increase is 1.5 to 4.5 times the global mean (Holland and Bitz, 2003). Warming-induced changes in carbon cycling are expected to exert large feedbacks on the global climate system (Davidson and Janssens, 2006; Christensen and Christensen, 2007; Oechel et al., 2000).

Warming is expected to accelerate soil C loss by increasing soil respiration but increasing nutrient mineralization, thereby stimulating plant net primary production (NPP) (Mack et al., 2004). Thus, the variation of climate may switch the role of the arctic system between a C sink and a source if soil C loss overtakes NPP (Davidson et al., 2000; Esteban G. Jobbágy and Jackson, 2000). Process-based biogeochemical models such as the Terrestrial Ecosystem Model (TEM; Hayes et al., 2014; Raich and Schlesinger, 1992; McGuire et al., 1992; Zhuang et al., 2001, 2002, 2003, 2010, 2013), Biome-BGC (Running and Coughlan, 1988), CASA (Potter et al., 1993), CENTURY (Parton et al., 1994), and the Biosphere Energy Transfer Hydrology scheme (BETHY) (Knorr, 2000) have been widely used to quantify the response of carbon dynamics to climatic changes (Todd-Brown et al., 2011). An ensemble of process-based model simulations suggests that arctic ecosystems acted as a sink of atmospheric CO₂ in recent decades (McGuire et al., 2012; Schimel, 2013). However, the response of this sink to increasing levels of atmospheric CO₂ and climate change is still uncertain (Todd-Brown et al., 2013). The IPCC fifth report also shows that land carbon storage is the largest source of uncertainty in the global carbon budget quantification (Ciais et al., 2013).

Much of the uncertainty is also due to the relatively lower levels of representation of ecosystem processes that determine the exchanges of water, energy, and C between land ecosystems and the atmosphere (Wieder et al., 2013) and ignorance of some key biogeochemical mechanisms (Schmidt et al., 2011). For example, heterotrophic respiration (R_H) is the primary loss pathway for soil organic carbon (Hanson et al., 2000; Bond-Lamberty and Thomson, 2010), and it generally increases with increasing temperature (Davidson and Janssens, 2006) and moisture levels in well-drained soils (Cook and Orchard, 2008). Moreover, this process is closely related to soil nitrogen mineralization that determines soil N availability and affects gross primary production (Hao et al., 2015). To date, most models treated soil decomposition as a first-order decay process, i.e., CO₂ respiration is directly proportional to soil organic carbon. However, it is not clear if these models are robust under changing environmental conditions (Lawrence et al., 2011; Schimel and Weintraub, 2003; Barichivich et al., 2013) since they often ignored the effects of changes in biomass and composition of decomposers, while recent empirical studies have shown that microbial abundance and community play a significant role in soil carbon decomposition (Allison and Martiny, 2008). The control that microbial activity and enzymatic kinetics imposed on soil respiration suggests the need for an explicit representation of microbial physiology, and enzymatic activity, in addition to the direct effects of soil temperature and soil moisture on heterotrophic respiration (Schimel and Weintraub, 2003). Recent mechanistically based models which explicitly incorporated the microbial dynamics and enzyme kinetics that catalyze soil C decomposition have produced notably different results and a closer match to contemporary observations (Wieder et al., 2013; Allison et al., 2010) indicating the need to incorporate these microbial mechanisms into large-scale earth system models to quantify carbon dynamics under future climatic conditions (Wieder et al., 2013; Allison et al., 2010).

This study advanced a microbe-based biogeochemistry model (MIC-TEM) based on the extant TEM (Raich and Schlesinger, 1992; McGuire et al., 1992; Zhuang et al., 2001, 2002, 2003, 2010, 2013; Hao et al., 2015). In MIC-TEM, the heterotrophic respiration is not only a function of soil temperature, soil organic matter (SOM), and soil moisture but also considers the effects of dynamics of microbial biomass and enzyme kinetics (Allison et al., 2010). The verified MIC-TEM was used to quantify the regional carbon dynamics in northern high latitudes (north of 45^∘ N) during the 20th and 21st centuries.

2 Methods

2.1 Overview

Below we first briefly describe how we advanced the MIC-TEM by modifying the soil respiration process in TEM (Zhuang et al., 2003) to better represent carbon dynamics in terrestrial ecosystems. Second, we describe how we parameterized and verified the new model using observed net ecosystem exchange data at representative sites and how simulated NPP was evaluated with Moderate Resolution Imaging Spectroradiometer (MODIS) data to demonstrate the reliability of a new model at regional scales. Third, we present how we applied the model to the northern high latitudes for the 20th and 21st centuries. Finally, we introduce how we conducted the sensitivity analysis on initial soil carbon input, using gridded observation-based soil carbon data of three soil depths during the 21st century.

2.2 Model description

TEM is a highly aggregated, large-scale biogeochemical model that estimates the dynamics of carbon and nitrogen fluxes and pool sizes of plants and soils using spatially explicit information on climate, elevation, soils, and vegetation (McGuire et al., 1992; Zhuang et al., 2003, 2010; Melillo et al., 1993). To explicitly consider the effects of microbial dynamics and enzyme kinetics on large-scale carbon dynamics of northern terrestrial ecosystems, we developed MIC-TEM by coupling version 5.0 of TEM (Zhuang et al., 2003, 2010) with a microbial-enzyme module (Hao et al., 2015; Allison et al., 2010). Our modification of the TEM improved the representation of the heterotrophic respiration (R_H) from a first-order structure to a more detailed structure (Fig. S1 in the Supplement).

In TEM, heterotrophic respiration R_H is calculated as a function of soil organic carbon (SOC), temperature sensitivity of heterotrophic soil respiration (Q₁₀), soil moisture (f(MOIST)), and the gram-specific decomposition constant K_d:

$$\begin{array}{}\text{(1)}& {\displaystyle }{R}_{\mathrm{H}}=K_{\mathrm{d}}\cdot \mathrm{SOC}\cdot Q_{\mathrm{10}}^{\frac{\mathrm{temp}.}{\mathrm{10}}}\cdot f\left(\mathrm{MOIST}\right),\end{array}$$

where temp is soil temperature in the top 20 cm (unit: ^∘C). CO₂ production from the SOC pool is directly proportional to the pool size, and the activity of decomposers only depends on the built-in relationships with soil temperature and moisture (Todd-Brown et al., 2012). Therefore, the changes in microbial community composition or adaption of microbial physiology to new conditions were not represented in TEM. However, current studies indicate that soil C decomposition depends on the activity of biological communities dominated by microbes (Schimel and Weintraub, 2003), implying that the biomass and composition of the decomposer community cannot be ignored (Todd-Brown et al., 2012).

We thus revised the first-order soil C structure in TEM to a second-order structure considering microbial dynamics and enzyme kinetics according to Allison et al. (2010). In MIC-TEM, heterotrophic respiration (R_H) is calculated as

$$\begin{array}{}\text{(2)}& {\displaystyle }{R}_{\mathrm{H}}=\mathrm{ASSIM}\cdot (\mathrm{1}-\mathrm{CUE}),\end{array}$$

where ASSIM and CUE represent microbial assimilation and carbon use efficiency, respectively. ASSIM is modeled with a Michaelis–Menten function:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{ASSIM}={\mathrm{Vmax}}_{\mathrm{uptake}}\cdot \mathrm{MIC}\cdot {\displaystyle \frac{\mathrm{DOC}}{{\mathrm{Km}}_{\mathrm{uptake}}+\mathrm{DOC}}},\end{array}$$

where DOC is dissolved organic carbon and Vmax_uptake is the maximum velocity of the reaction and calculated using the Arrhenius equation:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathrm{Vmax}}_{\mathrm{uptake}}={\mathrm{Vmax}}_{{\mathrm{uptake}}_{\mathrm{0}}}\cdot e^{-\frac{{\mathrm{Ea}}_{\mathrm{uptake}}}{\mathrm{R}\cdot \left(\mathrm{temp}+\mathrm{273}\right)}},\end{array}$$

where ${\mathrm{Vmax}}_{{\mathrm{uptake}}_{\mathrm{0}}}$ is the pre-exponential coefficient, Ea_uptake is the activation energy for the reaction (J mol⁻¹), R is the gas constant (8.314 J mol⁻¹ K⁻¹), and temp is the temperature in Celsius under the reaction occurs. Here we used soil temperature in the top 20 cm.

Moreover, Km_uptake is calculated as a linear function of temperature:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{Km}}_{\mathrm{uptake}}={\mathrm{Km}}_{{\mathrm{uptake}}_{\mathrm{slope}}}\cdot \mathrm{temp}+{\mathrm{Km}}_{{\mathrm{uptake}}_{\mathrm{0}}}.\end{array}$$

Microbial biomass MIC is modeled as

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{dMIC}}{\mathrm{d}t}}=\mathrm{ASSIM}\cdot \mathrm{CUE}-\mathrm{DEATH}-\mathrm{EPROD},\end{array}$$

where microbial biomass death (DEATH) and enzyme production (EPROD) are modeled as proportional to microbial biomass with rate constants r_death and r_EnzProd:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{DEATH}=r_{\mathrm{death}}\cdot \mathrm{MIC},\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{EPROD}=r_{\mathrm{EnzProd}}\cdot \mathrm{MIC},\end{array}$$

where r_death and r_EnzProd are the rate constants of microbial death and enzyme production, respectively.

DOC is part of soil organic carbon:

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{dDOC}}{\mathrm{d}t}}=& {\displaystyle }\mathrm{DEATH}\cdot \left(\mathrm{1}-\mathrm{MICtoSOC}\right)+\mathrm{DECAY}+\mathrm{ELOSS}\\ \text{(9)}& {\displaystyle }& {\displaystyle }-\mathrm{ASSIM},\end{array}$$

where MICtoSOC is carbon input ratio as dead microbial biomass to SOC, representing the fraction of microbial death that flows into SOC, and is set as a constant value according to Allison et al. (2010). SOC dynamics are modeled as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{dSOC}}{\mathrm{d}t}}=\mathrm{Litterfall}+\mathrm{DEATH}\cdot \mathrm{MICtoSOC}-\mathrm{DECAY},\end{array}$$

where litterfall is estimated as a function of vegetation carbon (Zhuang et al., 2010). The enzymatic decay of SOC is calculated as

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathrm{DECAY}=\mathrm{Vmax}\cdot \mathrm{ENZ}\cdot {\displaystyle \frac{\mathrm{SOC}}{\mathrm{Km}+\mathrm{SOC}}},\end{array}$$

where Vmax is the maximum velocity of the reaction and calculated using the Arrhenius equation:

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathrm{Vmax}={\mathrm{Vmax}}_{\mathrm{0}}\cdot e^{-\frac{\mathrm{Ea}}{\mathrm{R}\cdot \left(\mathrm{temp}+\mathrm{273}\right)}}.\end{array}$$

The parameters Km and carbon use efficiency (CUE) are temperature sensitive and calculated as a linear function of temperature between 0 and 50 ^∘C:

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{Km}={\mathrm{Km}}_{\mathrm{slope}}\cdot \mathrm{temp}+{\mathrm{Km}}_{\mathrm{0}}\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{CUE}={\mathrm{CUE}}_{\mathrm{slope}}\cdot \mathrm{temp}+{\mathrm{CUE}}_{\mathrm{0}},\end{array}$$

where CUEslope and CUE₀ are parameters for calculating CUE and Km_slope and Km₀ are parameters for calculating Km. The values of CUE_slope, CUE₀, Km_slope, and Km₀ were derived from Allison et al. (2010).

ELOSS is also a first-order process, representing the loss of enzymes:

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{ELOSS}=r_{\mathrm{enzloss}}\cdot \mathrm{ENZ},\end{array}$$

where r_enzloss is the rate constant of enzyme loss. The enzyme pool (ENZ) is modeled as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{dENZ}}{\mathrm{d}t}}=\mathrm{EPROD}-\mathrm{ELOSS}.\end{array}$$

Heterotrophic respiration (R_H) is an indispensable component of soil respiration (Bond-Lamberty and Thomson, 2010) and is closely coupled with soil nitrogen (N) mineralization, which determines soil N availability, affecting gross primary production (GPP).

2.3 Model parameterization and validation

Table 1 Parameters associated with more detailed microbial dynamics in MIC-TEM.

The variables and parameters of these microbial dynamics and their impacts on soil C decomposition were detailed in Allison et al. (2010) (Table 1). Here we parameterized MIC-TEM for representative ecosystem types in northern high latitudes based on monthly net ecosystem production (NEP, gCm⁻² mon⁻¹) measurements from the AmeriFlux network (Davidson et al., 2000) (Table S1 in the Supplement). The results for model parameterization are presented in Fig. S2. Another set of level-4 gap-filled NEP data were used for model validation at site level (Table S2). The site-level monthly climate data of air temperature (^∘C), precipitation (mm), and cloudiness (%) were used to drive the model. Gridded MODIS NPP data from 2001 to 2010 were used to evaluate regional NPP simulations. The MODIS NPP data were developed by the MOD17 MODIS project. The product name is Net Primary Production Yearly L4 Global 1 km. The critical parameter used in the MOD17 algorithm is conversion efficiency parameter ε. More information about the MODIS NPP product can be found at https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD17A2_M_PSN (last access: January 2017).

In TEM, NPP is calculated as

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathrm{NPP}=\mathrm{GPP}-R_{\mathrm{A}},\end{array}$$

where GPP is gross primary production and R_A is autotrophic respiration. GPP is defined as

$$\begin{array}{ll}{\displaystyle }\mathrm{GPP}=& {\displaystyle }{C}_{max}\cdot f\left(\mathrm{PAR}\right)\cdot f\left(\mathrm{phenology}\right)\cdot f\left(\mathrm{foliage}\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }\cdot f\left(T\right)\cdot f\left({\mathrm{CO}}_{\mathrm{2}}\right)\cdot f\left(\mathrm{NA}\right)\cdot f\left(\mathrm{FT}\right),\end{array}$$

where C_max is the maximum rate of carbon assimilation, PAR is photosynthetically active radiation, and f(phenology) represents the effects of leaf area (Raich and Schlesinger, 1992). The function f(foliage) represents the ratio of canopy leaf biomass relative to maximum leaf biomass (Zhuang et al., 2002). T is monthly air temperature, and f(CO₂) represents the effects of elevated atmospheric CO₂ (McGuire et al., 1992; Pan et al., 1998). The function f(NA) models the limiting effects of plant nitrogen status on GPP (McGuire et al., 1992; Pan et al., 1998). The function f(FT) represents the effects of freeze–thaw (Zhuang et al., 2003). For detailed GPP and R_A calculations, see Zhuang et al. (2003).

The parameterization was conducted with a global optimization algorithm SCE-UA (Shuffled Complex Evolution) (Duan et al., 1994) to minimize the difference between the monthly simulated and measured net ecosystem exchange (NEE) at these sites (Fig. S2). The cost function of the minimization is

$$\begin{array}{}\text{(19)}& {\displaystyle }\mathrm{Obj}=\sum _{i=\mathrm{1}}^k({\mathrm{NEP}}_{\mathrm{obs},i}-{\mathrm{NEP}}_{\mathrm{sim},i}{)}^{\mathrm{2}},\end{array}$$

where NEP_obs,i and NEP_sim,i are the observed and simulated NEP, respectively. k is the number of data pairs for comparison. Other parameters used in MIC-TEM were default values from TEM 5.0 (Zhuang et al., 2003, 2010). The optimized parameters were used for model validation and regional extrapolations.

2.4 Regional simulations

Two sets of regional simulations for the 20th century using MIC-TEM and TEM at a spatial resolution of 0.5^∘ latitude × 0.5^∘ longitude were conducted. Gridded forcing data of monthly air temperature, precipitation, and cloudiness were used, along with other ancillary inputs including historical atmospheric CO₂ concentrations, soil texture, elevation, and potential natural vegetation. Climatic inputs vary over time and space, whereas soil texture, elevation, and land cover data are assumed to remain unchanged throughout the 20th century and only vary spatially. The transient climate data during the 20th century were obtained from the Climatic Research Unit (CRU TS3.1) of the University of East Anglia (Harris et al., 2014). The spatially explicit data include potential natural vegetation (Melillo et al., 1993), soil texture (Zhuang et al., 2003), and elevation (Zhuang et al., 2015).

Similarly, two sets of simulations were conducted driven with two contrasting climate change scenarios (RCP 2.6 and RCP 8.5) over the 21st century. The future climate change scenarios were derived from the HadGEM2-ES model, which is a member of the CMIP5 project (https://esgf-node.llnl.gov/search/cmip5/, January 2017). The future atmospheric CO₂ concentrations and climate forcing from each of the two climate change scenarios were used. The simulated NPP, R_H, and NEP by both models (TEM 5.0 and MIC-TEM) were analyzed. The positive NEP represents a CO₂ sink from the atmosphere to terrestrial ecosystems, while a negative value represents a source of CO₂ from terrestrial ecosystems to the atmosphere.

Moreover, in order to test the parameter uncertainty in our model, we conducted the regional simulations with 50 sets of parameters for both historical and future studies. The 50 sets of parameters were obtained according to the method in Tang and Zhuang (2008). The upper and lower bounds of the regional estimations were generated based on these simulations.

2.5 Sensitivity to initial soil carbon input

Future carbon dynamics can be affected by varying initial soil carbon amount. In the standard simulation of TEM, the initial soil carbon amount for transient simulations was obtained from equilibrium and spin-up periods directly for each grid cell in the region. To test the sensitivity to the initial soil carbon amount in transient simulations for the 21st century, we used empirical soil organic carbon data extracted from the Northern Circumpolar Soil Carbon Database (NCSCD) (Tarnocai et al., 2009), as the initial soil carbon amount. The ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$ soil carbon data products for three different depths of 30, 100, and 300 cm were used. The sensitivity test was conducted for transient simulations under the RCP 2.6 and RCP 8.5 scenarios. To avoid the instability of the C : N ratio caused by replacing the initial soil carbon pool with observed data at the beginning of transient period, initial soil nitrogen values were also generated based on the soil carbon data and corresponding C : N ratio map for transient simulations (Zhuang et al., 2003; Raich and Schlesinger, 1992).

3 Results

3.1 Model verification at site and regional levels

With the optimized parameters, MIC-TEM reproduces the carbon dynamics well for alpine tundra, boreal forest, temperate coniferous forest, temperate deciduous forest, grasslands, and wet tundra with R² ranging from 0.70 for Ivotuk to 0.94 for the Bartlett Experimental Forest (Fig. S3, Table S3). In general, the model performs better for forest ecosystems than for tundra ecosystems. The temporal NPP from 2001 to 2010 simulated by MIC-TEM and TEM was compared with MODIS NPP data (Fig. S4). Pearson correlation coefficients are 0.52 (MIC-TEM and MODIS) and 0.34 (TEM and MODIS). NPP simulated by MIC-TEM showed higher spatial correlation coefficients with MODIS data than TEM (Fig. S5). By considering more detailed microbial activities, the heterotrophic respiration is more adequately simulated using the MIC-TEM. The simulated differences in soil decomposition result in different levels of soil-available nitrogen, which influences the nitrogen uptake by plants, the rate of photosynthesis, and NPP. The spatial correlation coefficient between NPP simulated by MIC-TEM and MODIS is close to 1 in most study areas, suggesting the reliability of MIC-TEM at the regional scale.

3.2 Regional carbon dynamics during the 20th century

Figure 1 Simulated annual net primary production (NPP, a), heterotrophic respiration (R_H, b), and net ecosystem production (NEP, c) by MIC-TEM with an ensemble of parameters.

The equifinality of the parameters in MIC-TEM was considered in our ensemble regional simulations to measure the parameter uncertainty (Tang and Zhuang, 2008). Here and below, the ensemble means and the inter-simulation standard deviations are shown to measure the uncertainty unless otherwise specified. These ensemble simulations indicated that the northern high latitudes go from acting as a carbon source of 38.9 Pg C to being a carbon sink of 190.8 Pg C by different ensemble members, with a mean of 64.2±21.4 Pg at the end of the 20th century, while the simulation with the optimized parameters estimates a regional carbon sink of 77.6 Pg with an interannual standard deviation of 0.21 Pg C yr⁻¹ during the 20th century (Fig. 1). Simulated regional NEP with optimized parameters using TEM and MIC-TEM showed an increasing trend throughout the 20th century except for a slight decrease during the 1960s (Fig. 2). The spatial distributions of NEP simulated by MIC-TEM for different periods in the 20th century also show an increasing trend (Fig. 3). Positive values of NEP represent sinks of CO₂ into terrestrial ecosystems, while negative values represent sources of CO₂ to the atmosphere. From 1900 onwards, both models estimated a regional carbon sink during the 20th century. With optimized parameters, TEM estimated higher NPP and R_H at 0.6 and 0.3 Pg C yr⁻¹ than MIC-TEM, respectively, at the end of the 20th century (Fig. 2). MIC-TEM estimated a carbon sink increase from 0.64 to 0.83 Pg C yr⁻¹ during the century while the estimated increase by TEM was much higher (0.28 Pg C yr⁻¹) (Fig. 2). At the end of the century, MIC-TEM-estimated NEP reached 1.0 Pg C yr⁻¹ in comparison with TEM estimates of 0.3 Pg C yr⁻¹. TEM-estimated NPP and R_H are 0.5 and 0.3 Pg C yr⁻¹ higher, respectively. As a result, TEM estimated that the region accumulated 11.4 Pg more carbon than MIC-TEM. Boreal forests are a major carbon sink at 0.55 and 0.63 Pg C yr⁻¹ estimated by MIC-TEM and TEM, respectively. Alpine tundra contributes the least sink. Overall, TEM overestimated the sink by 12.5 % in comparison to MIC-TEM for forest ecosystems and 16.7 % for grasslands. For wet tundra and alpine tundra, TEM overestimated about 20 % and 33 %, respectively, in comparison with MIC-TEM (Table 2).

Figure 2 Simulated annual net primary production (NPP, a), heterotrophic respiration (R_H, b), and net ecosystem production (NEP, c) by MIC-TEM and TEM.

Table 2 Partitioning of average annual net ecosystem production (as Pg C per year) for six vegetation types during the 20th century.

3.3 Regional carbon dynamics during the 21st century

Simulated regional annual NPP and R_H increase under the RCP 8.5 scenario with both models (Fig. 4). With optimized parameters, MIC-TEM-estimated NPP increases from 9.2 in the 2000s to 13.2 Pg C yr⁻¹ in the 2090s, while TEM-predicted NPP is 2.0 Pg C yr⁻¹ higher in the 2000s and 0.3 Pg C yr⁻¹ higher in the 2090s (Fig. 4). Similarly, TEM also overestimated R_H by 1.7 Pg C yr⁻¹ in the 2000s and by 0.25 Pg C yr⁻¹ in the 2090s, respectively (Fig. 4). As a result, the regional sink increases from 0.53 Pg C yr⁻¹ in the 2000s to 1.4 Pg C yr⁻¹ in the 2070s and then decreases to 1.1 Pg C yr⁻¹ in the 2090s as estimated by MIC-TEM (Fig. 4). Given the uncertainty in parameters, MIC-TEM predicted that the region acts as a carbon sink of 48.7 to 140.7 Pg, with a mean of 71.7±26.6 Pg at the end of 21st century, while the simulation with optimized parameters estimates a regional carbon source of 79.5 Pg with an interannual standard deviation of 0.37 Pg C yr⁻¹ during the 21st century (Fig. 4). TEM predicted a similar trend for NEP, which overestimated the carbon sink with a magnitude of 19.2 Pg compared to the simulation by MIC-TEM with optimized parameters. Under the RCP 2.6 scenario (Fig. 4), the increase in NPP and R_H is smaller from 2000 to 2100 compared to the simulation under the RCP 8.5. MIC-TEM predicted that NPP increases from 9.1 to 10.9 Pg C yr⁻¹, and TEM estimated it to be 1.6 Pg C yr⁻¹ higher at the beginning and 0.9 Pg C yr⁻¹ higher in the end of the 21st century (Fig. 4). Consequently, MIC-TEM-predicted NEP fluctuates between sinks and sources during the century, being neutral before 2070 and a source of between −0.2 and −0.3 Pg C yr⁻¹ after the 2070s. As a result, the region acts as a carbon source of 1.6 Pg C with an interannual standard deviation of 0.24 Pg C yr⁻¹ estimated with MIC-TEM and a sink of 27.6 Pg C with an interannual standard deviation of 0.2 Pg C yr⁻¹ estimated with TEM during the century (Fig. 4). When considering the uncertainty source of parameters, MIC-TEM predicted that the region goes from acting as a carbon source of 64.8 Pg C to being a carbon sink of 58.6 Pg C during the century with a mean of $-\mathrm{3.3}\pm \mathrm{20.3}$ Pg at the end of 21st century (Fig. 4).

Figure 3 Spatial distribution of NEP simulated by MIC-TEM for the periods (a) 1900–1930, (b) 1931–1960, (c) 1961–1990, and (d) 1991–2000. Positive values of NEP represent sinks of CO₂ into terrestrial ecosystems, while negative values represent sources of CO₂ to the atmosphere.

Figure 4 Predicted changes in carbon fluxes: (a, b) NPP, (c, d) R_H, and (e, f) NEP for all land areas north of 45^∘ N in response to transient climate change under (a, c, e) the RCP 8.5 scenario and (b, d, f) the RCP 2.6 scenario with MIC-TEM and TEM 5.0. The decadal running mean is applied. The blue area represents the upper and lower bounds of simulations.

3.4 Model sensitivity to initial soil carbon

Under the RCP 2.6, without replacing the initial soil carbon with inventory-based estimates (Tarnocai et al., 2009) in model simulations, TEM estimated that the regional SOC is 604.2 Pg C and accumulates an additional 12.1 Pg C during the 21st century. When using estimated soil carbon (Tarnocai et al., 2009) within depths of 30, 100, and 300 cm as initial pools in simulations, TEM predicted that regional SOC is 429.5, 689.3, and 1003.4 Pg C in 2000 and increases by 9.9, 16.0, and 22.8 Pg C at the end of the 21st century, and the regional cumulative carbon sink is 20.4, 34.0, and 48.1 Pg C, respectively, during the century. In contrast, using the same inventory-based SOC estimates, MIC-TEM projected that the region goes from acting as a cumulative carbon sink to being a source at 0.7, 2.2, and 3.0 Pg C. Under the RCP 8.5, both models predicted that the region acts as a carbon sink, regardless of the magnitudes of initial soil carbon pools used, with a TEM-projected sink of 71.7, 120, and 155.6 Pg C and a much smaller cumulative sink of 65.4, 88.6, and 109.8 Pg C estimated with MIC-TEM (Table 3).

Table 3 Increasing of SOC, vegetation carbon (VGC), soil organic nitrogen (SON), vegetation nitrogen (VGN) from 1900 to 2000, and total carbon storage during the 21st century predicted by two models with observed soil carbon data of three different depths under (a) RCP 2.6 and (b) RCP 8.5.

4 Discussion

During the last few decades, a greening accompanying warming and rising atmospheric CO₂ in the northern high latitudes (>45^∘ N) has been documented (McGuire et al., 1995; McGuire and Hobbie, 1997; Chapin and Starfield, 1997; Stow et al., 2004; Callaghan et al., 2005; Tape et al., 2006). The large stocks of carbon contained in the region (Tarnocai et al., 2009) are particularly vulnerable to climate change (Schuur et al., 2008; McGuire et al., 2009). To date, the degree to which the ecosystems may serve as a source or a sink of C in the future is still uncertain (McGuire et al., 2009; Wieder et al., 2013). Therefore, accurate models are essential for predicting carbon–climate feedbacks in the future (Todd-Brown et al., 2013). Our regional simulations indicate that the region is currently a carbon sink, which is consistent with many previous studies (White et al., 2000; Houghton, 2007), and this sink will grow under the RCP 8.5 scenario but shift to a carbon source under the RCP 2.6 scenario by 2100. MIC-TEM shows a higher correlation between NPP and soil temperature (R=0.91) than TEM (R=0.82), suggesting that MIC-TEM is more sensitive to environmental changes (Table S4).

Our regional estimates of carbon fluxes by MIC-TEM are within the uncertainty range from other existing studies. For instance, Zhuang et al. (2003) estimated the region to be a sink of 0.9 Pg C yr⁻¹ in extratropical ecosystems for the 1990s, which is similar to our estimation of 0.83 Pg C yr⁻¹ by MIC-TEM. White et al. (2000) estimated that, during the 1990s, regional NEP above 50^∘ N is 0.46 Pg C yr⁻¹ while Qian et al. (2010) estimated that NEP increased from 0 to 0.3 Pg C yr⁻¹ for the high-latitude region above 60^∘ N during last century and reached 0.25 Pg C yr⁻¹ during the 1990s. White et al. (2000) predicted that, from 1850 to 2100, the region accumulated 134 Pg C in terrestrial ecosystems, in comparison with our estimates of 77.6 Pg C with MIC-TEM and 89 Pg C with TEM. Our projection of a weakening sink during the second half of the 21st century is consistent with previous model studies (Schaphoff et al., 2013). Our predicted trend of NEP is very similar to the finding of White et al. (2000), indicating that NEP increases from 0.46 Pg C yr⁻¹ in the 2000s and reaches 1.5 Pg C yr⁻¹ in the 2070s, then decreases to 0.6 Pg C yr⁻¹ in the 2090s.

The MIC-TEM-simulated NEP generally agrees with the observations. However, model simulations still deviate from the observed data, especially for tundra ecosystems. The deviation may be due to the uncertainty or errors in the observed data, which do not constrain the model parameters well. Uncertain driving data such as temperature and precipitation are also a source of uncertainty for transient simulations. In addition, we assumed that vegetation will not change during the transient simulation. However, over the past few decades in the northern high latitudes, temperature increases have led to vegetation changes (Hansen et al., 2006), including latitudinal tree line advance (Lloyd, 2005) and increasing shrub density (Sturm et al., 2001). Vegetation can shift from one type to another because of competition for light, N, and water (White et al., 2000). For example, needleleaf trees tend to replace tundra gradually in response to warming. In some areas, forests even moved several hundreds of kilometers within 100 years (Gear and Huntley, 1991). The vegetation changes will affect carbon cycling in these ecosystems. In addition, we have not yet considered the effects of management of agricultural lands (Cole et al., 1997), but Zhuang et al. (2003) showed that the changes in agricultural land use in northern high latitudes have been small.

The largest limitation on this study is that we have not explicitly considered the fire effects. Warming in the northern high latitudes could favor fire in its frequency, intensity, seasonality, and extent (Kasischke and Turetsky, 2006; Johnstone and Kasischke, 2005; Soja et al., 2007; Randerson et al., 2006; Bond-Lamberty et al., 2007). Fire has profound effects on northern forest ecosystems, altering the N cycle and water and energy exchanges between the atmosphere and ecosystems. An increase in wildfires will destroy most of aboveground biomass and consume organic soils, resulting in less carbon uptake by vegetation (Harden et al., 2000), leading to a net release of carbon in the short term. However, a suite of biophysical mechanisms of ecosystems, including post-fire increase in the surface albedo and rates of biomass accumulation, may, in turn, exert a negative feedback on climate warming (Amiro et al., 2006; Goetz et al., 2007), further influencing the carbon exchanges between ecosystems and the atmosphere.

Moreover, carbon uptake in land ecosystems depends on new plant growth, which connects closely with the availability of nutrients such as mineral nitrogen. Recent studies have shown that when soil nitrogen is in short supply, most terrestrial plants form symbiotic relationships with fungi; hyphae provide nitrogen to plants, and in return, plants provide sugar to fungi (Hobbie and Hobbie, 2008, 2006; Schimel and Hättenschwiler, 2007). This symbiotic relationship has not been considered in our current modeling, which may lead to a large uncertainty in our quantification of carbon and nitrogen dynamics.

A shift in microbial community structure was not considered in our model, which could affect the temperature sensitivity of heterotrophic respiration (Stone et al., 2012). The Michaelis–Menten constant (Km) could also adapt to climate warming, and it may increase more significantly with increasing temperature in cold-adapted enzymes than in warm-adapted enzymes (German et al., 2012; Somero, 2004; Dong and Somero, 2009). Carbon use efficiency (CUE) is also a controversial parameter in our model. Empirical studies in soils suggest that microbial CUE declines by at least 0.009 ^∘C⁻¹ (Steinweg et al., 2008), while other studies find that CUE is invariant with temperature (López-Urrutia and Morán, 2007). Another key microbial trait lacking in our modeling is microbial dormancy (He et al., 2015). Dormancy is a common, bet-hedging strategy used by microorganisms when environmental conditions limit their growth and reproduction (Lennon and Jones, 2011). Microorganisms in dormancy are not able to drive biogeochemical processes such as soil CO₂ production, and, therefore, only active microorganisms should be involved in utilizing substrates in soils (Blagodatskaya and Kuzyakov, 2013). Many studies have indicated that soil respiration responses to environmental conditions are more closely associated with the active portion of microbial biomass than total microbial biomass (Hagerty et al., 2014; Schimel and Schaeffer, 2012; Steinweg et al., 2013). Thus, the ignorance of microbial dormancy could fail to distinguish between microbes with different physiological states, introducing uncertainties into our carbon estimation.

5 Conclusions

This study used a more detailed microbial biogeochemistry model to investigate the carbon dynamics in the region for the past and this century. Regional simulations using MIC-TEM indicated that, over the 20th century, the region is a sink of 77.6 Pg C. This sink could reach 79.5 Pg C under the RCP 8.5 scenario or shift to a carbon source of 1.6 Pg under the RCP 2.6 scenario during the 21st century. On the other hand, traditional TEM overestimated the carbon sink under the RCP 8.5 scenario with a magnitude of 19.2 Pg than MIC-TEM and predicted that this region would act as a carbon sink with a magnitude of 27.6 Pg under the RCP 2.6 scenario during the 21st century. Using recent soil carbon stock data as initial soil carbon in model simulations, the region was estimated to shift from a carbon sink to a source, with total carbon release at 0.7–3 Pg by 2100 depending on initial soil carbon pools at different soil depths under the RCP 2.6 scenario. In contrast, the region acts as a carbon sink at 55.4–99.8 Pg in the 21st century under the RCP 8.5 scenario. Without considering more detailed microbial processes, models estimated that the region acts as a carbon sink under both scenarios. Under the RCP 2.6 scenario, the cumulative sink ranges from 9.9 to 22.8 Pg C. Under the RCP 8.5 scenario, the cumulative sink is even larger at 71.7–155.6 Pg C. This study indicated that more detailed microbial physiology-based biogeochemistry models estimate carbon dynamics very differently from using a relatively simple microbial decomposition-based model. The comparison with satellite products or other estimates for the 20th century suggests that the more detailed microbial decomposition should be considered to adequately quantify C dynamics in northern high latitudes.