Microbial dormancy and its impacts on northern temperate and boreal terrestrial ecosystem carbon budget

A large amount of soil carbon in northern temperate and boreal regions could be emitted as greenhouse gases in a warming future. However, lacking detailed microbial processes such as microbial dormancy in current biogeochemistry models might have biased the quantification of the regional carbon dynamics. Here the effect of microbial dormancy was incorporated into a biogeochemistry model to improve the quantification for the last century and this century. Compared with the previous model without considering the microbial dormancy, the new model estimated the regional soils stored 75.9 Pg more C in the terrestrial ecosystems during the last century and will store 50.4 and 125.2 Pg more C under the RCP8.5 and RCP2.6 scenarios, respectively, in this century. This study highlights the importance of the representation of microbial dormancy in earth system models to adequately quantify the carbon dynamics in the northern temperate and boreal natural terrestrial ecosystems.

1 Introduction

The land ecosystems in northern temperate and boreal regions (>45^∘ N) occupy 22 % of the global surface and store over 40 % of the global soil organic carbon (SOC) (McGuire and Hobbie, 1997; Melillo et al., 1993; Tarnocai et al., 2009; Hugelius et al., 2014). During the past decades, a greening accompanying a warming in the region has been documented (Zhou et al., 2001; Lloyd et al., 2002; Stow et al., 2004; Callaghan et al., 2005; Tape et al., 2006). The regional carbon dynamics are expected to loom large in the global carbon cycle and exert large feedbacks to the global climate system (McGuire et al., 2009; Davidson and Janssens, 2006; Bond-Lamberty and Thomson, 2010).

To date, numerous ecosystem models have been developed to project the feedbacks between terrestrial ecosystem carbon cycling and climate (Raich et al., 1991; Zhuang et al., 2001, 2002, 2015; Parton et al., 1993; Knorr et al., 2005; Running and Coughlan, 1988), but they can bias their quantifications due to missing detailed microbial mechanisms in these models (Schmidt et al., 2011; Todd-Brown et al., 2013; Conant et al., 2011; Treseder et al., 2011). Microorganisms play a central role in decomposition of litter and soil organic carbon, which further governs the global carbon cycling and climate change (Xu et al., 2014; Treseder et al., 2011; Wang et al., 2015). An emerging field of research has begun to incorporate microbial ecology into existing process-based models to represent decomposition in ways that include important microbial processes that were previously ignored (Zha and Zhuang, 2018; Schimel and Weintraub, 2003; Allison et al., 2010; German et al., 2012). These microbe-based models tend to better reproduce field and satellite observations than traditional ones that treat soil decomposition as a first-order decay process without considering microbial activities (Treseder et al., 2011; Wieder et al., 2013; Todd-Brown et al., 2011; Lawrence et al., 2009; Moorhead and Sinsabaugh, 2006). However, some vital microbial traits such as microbial dormancy and community shifts are still rarely explicitly considered in large-scale ecosystem models (Wieder et al., 2015), and this may introduce notable uncertainties (Graham et al., 2014, 2016; Wang et al., 2015; Bouskill et al., 2012; Kaiser et al., 2014).

Dormancy is broadly recognized as a strategy for microorganisms to cope with periodical environmental stresses (Harder and Dijkhuizen, 1983). When environmental conditions are unfavorable for growth, microbes switch to a dormant state, which is a reversible state of low to zero metabolic activity (Stolpovsky et al., 2011; Lennon and Jones, 2011). In this state, biogeochemical processes such as soil decomposition are slow (Blagodatskaya and Kuzyakov, 2013). At any given time, there is only a fraction of, likely below 50 %, metabolically active microbes in natural soils (Wang et al., 2015; Stolpovsky et al., 2011). Soil decomposition and nutrient cycling mainly depend on these active microbes because only active ones can consume organic matter and replicate themselves (Wang et al., 2015; Blagodatskaya et al., 2014). To date, most existing biogeochemistry models use total rather than active microbial biomass as an indicator of microbial activities (Wieder et al., 2015), which could bias the estimates of soil decomposition and ecosystem carbon budget (Hagerty et al., 2014; He et al., 2015). In particular, the northern temperate and boreal terrestrial ecosystems are nitrogen-limited; neglecting microbial dormancy will lead to incorrect estimates of nitrogen availability through soil decomposition, failing to capture nitrogen feedbacks to carbon dynamics (Wang et al., 2015; Stolpovsky et al., 2011; Thullner et al., 2005). Furthermore, these ecosystems have experienced a marked seasonality of active and dormant microbial cycles and the above-global-average warming, which might have increased the proportion of active microbes in soils (He et al., 2015). Thus, incorporating dormancy effects will improve model realism to provide a better projection of the northern temperate and boreal terrestrial ecosystem carbon dynamics.

This study incorporated the effects of the microbial dormancy trait into an extant process-based biogeochemistry model (MIC-TEM) (Zha and Zhuang, 2018; He et al., 2015). The dormant and active microbial physiology has been considered explicitly in the new version of the model (MIC-TEM-dormancy). The revised model was parameterized, validated, and then applied to evaluate the carbon dynamics during the last century and this century in the northern temperate and boreal terrestrial ecosystems (north 45^∘ N above). By comparing the results of MIC-TEM-dormancy and MIC-TEM, we can show that incorporating microbial dormancy may produce a much different prediction in the historical and future carbon budget.

2 Methods

2.1 Overview

Due to the importance of microbial dormancy, some recent work has been done to consider the metabolic activation and deactivation of microbes in soil and its effects on soil carbon (C) dynamics and climate feedbacks. For example, Wang et al. (2015) have incorporated transformation processes between active and dormant states to develop two versions of MEND, that is, MEND with and without dormancy. The two versions of the model have been applied to quantify the carbon decomposition in laboratory incubations of four soils. Salazar et al. (2018) have also taken microbially dormancy into account to compare their predictions of microbial biomass and soil heterotrophic respiration (R_H) under simulated cycles of stressful (dryness) and favorable (wet pulses) conditions. Our study extends those modeling studies to the northern temperate and boreal terrestrial ecosystems by developing a more detailed biogeochemistry model considering the dormancy impacts. Below, we first describe how we developed the new model (MIC-TEM-dormancy) by incorporating the microbial dormancy trait into an existing microbe-based biogeochemistry model (MIC-TEM). Second, we discuss how parameterization and validation of the MIC-TEM-dormancy model were conducted using observed net ecosystem exchange data and heterotrophic respiration data at representative sites. Third, we presented how the model was applied to natural ecosystems in the region (above 45^∘ N) for the 20th and 21st centuries and discussed the dormancy effects on their regional carbon budget.

2.2 Model description

A non-dormancy version of the biogeochemistry model (MIC-TEM) has been developed by incorporating a microbial module (Allison et al., 2010) into an extant large-scale biogeochemical model (TEM) to explicitly (Zhuang et al., 2015) consider the effects of microbial dynamics and enzyme kinetics on carbon dynamics (Zha and Zhuang, 2018). Here we further advanced the MIC-TEM by incorporating algorithms that describe the effects of microbial dormancy dynamics based on He et al. (2015). Different from He et al. (2015), in which the microbial module was driven with existing data of carbon stocks and fluxes, our study incorporated the microbial module into an extant MIC-TEM that simulates carbon data dynamically. This coupling enables us to extrapolate our model to northern temperate and boreal terrestrial ecosystems, rather than only for temperate forest regions in He et al. (2015). In our new model (MIC-TEM-dormancy), the microbial biomass pool was divided into two fractions, including the dormant and active microbial biomass pools. The two microbial biomass pools and the reversible transition between them have been considered explicitly in the new model (Fig. 1) but were ignored in MIC-TEM.

Figure 1 Framework of the dormancy model: microbial biomass is split into two parts, active microbial biomass and dormant microbial biomass (shown in the green dashed circle). Maintenance respiration from these two parts and the CO₂ production through microbial assimilation contribute to heterotrophic respiration. The model was revised based on Zha and Zhuang (2018).

In previous MIC-TEM, heterotrophic respiration (R_H) is simply calculated as the product of ASSIM and CUE, which are microbial assimilation and carbon use efficiency, respectively. For detailed carbon dynamics in MIC-TEM, see Zha and Zhuang (2018).

Here we revised MIC-TEM by incorporating microbial dormancy dynamics according to He et al. (2015). In MIC-TEM-dormancy, the soil heterotrophic respiration R_H is revised to include three parts: the maintenance respiration from the active and dormant microorganisms and the CO₂ production through the process of microbial assimilation (He et al., 2015):

$$\begin{array}{}\text{(1)}& R_{\mathrm{H}}=m_{\mathrm{R}}{Q}_{\text{10 mic}}^{\frac{\text{temp-15}}{\mathrm{10}}}{B}_{\text{a}}+\mathit{\beta }{\mathrm{m}}_{\mathrm{R}}{Q}_{\text{10 mic}}^{\frac{\text{temp-15}}{\mathrm{10}}}{B}_{\mathrm{d}}+{\text{CO}}_{\mathrm{2}},\end{array}$$

where the first two terms are maintenance respiration from the active and dormant microorganisms. The last term is the CO₂ produced during the process of microbial assimilation.

For the first two terms, B_a and B_d represent the active and dormant microbial biomass pool, respectively. The parameter m_R denotes the specific maintenance rate in an active state (h⁻¹), and β is the ratio of dormant maintenance rate to active maintenance rate. Thus, βm_R denotes the maximum specific maintenance rate in a dormant state. Temperature sensitivity was expressed as the Q₁₀ function $\left(Q_{\mathrm{10}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}\right)$, where temp is soil temperature in the top 20 cm (degrees Celsius).

For the third term, the CO₂ produced through microbial assimilation is calculated as in He et al. (2015) and Allison et al. (2010):

$$\begin{array}{}\text{(2)}& {\mathrm{CO}}_{\mathrm{2}}=\text{ASSIM}\times (\mathrm{1}-Y_{\mathrm{g}}),\end{array}$$

where ASSIM represents the microbial assimilation, and the parameter Y_g represents carbon use efficiency. Microbial assimilation (ASSIM) is calculated as in He et al. (2015):

$$\begin{array}{}\text{(3)}& \mathrm{ASSIM}={\displaystyle \frac{\mathrm{1}}{Y_{\mathrm{g}}}}{\displaystyle \frac{\mathrm{\Phi }}{\mathit{\alpha }}}{m}_{\mathrm{R}}{Q}_{\mathrm{10}\mathrm{enz}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}{B}_{\mathrm{a}}{\left({\displaystyle \frac{{\mathrm{CN}}_{\mathrm{soil}}}{{\mathrm{CN}}_{\mathrm{mic}}}}\right)}^{\mathrm{0.6}}.\end{array}$$

Here parameter α is maintenance weight (h⁻¹); CN_soil and CN_mic denote the C:N ratios of soil and that of microbial biomass. In addition, Φ is the substrate saturation level and defined as in He et al. (2015) and Wang et al. (2014):

$$\begin{array}{}\text{(4)}& \mathrm{\Phi }={\displaystyle \frac{\mathrm{S}}{K_{\mathrm{s}}+\mathrm{S}}},\end{array}$$

where K_s is the half saturation constant for substrate uptake as indicated by the Michaelis–Menten kinetic, and S is soluble C substrates that are directly accessible for microbial assimilation (Wang et al., 2014). Here we quantified concentration of soluble C substrates that are directly accessible for microbial assimilation by using the conceptual framework from Davidson et al. (2012):

$$\begin{array}{}\text{(5)}& \mathrm{S}=\mathrm{soluble}\mathrm{C}\times D_{\mathrm{liq}}\times {\mathit{\theta }}^{\mathrm{3}}.\end{array}$$

The term “soluble C” denotes the state variable of the soluble carbon pool. D_liq is the diffusion coefficient of the substrate in the liquid phase and is formulated as

$$\begin{array}{}\text{(6)}& D_{\mathrm{liq}}=\mathrm{1}/(\mathrm{1}-\mathrm{BD}/\mathrm{PD}{)}^{\mathrm{3}},\end{array}$$

where BD is the bulk density and PD is the soil particle density. θ is the volumetric soil moisture.

Different from MIC-TEM, the transitions between active and dormant microbial biomass are included in MIC-TEM-dormancy.

$$\begin{array}{}\text{(7)}& {\displaystyle }{B}_{\mathrm{a}\to \mathrm{d}}=(\mathrm{1}-\mathrm{\Phi })m_{\mathrm{R}}{Q}_{\mathrm{10}\mathrm{mic}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}{B}_{\mathrm{a}},\text{(8)}& {\displaystyle }{B}_{\mathrm{d}\to \mathrm{a}}=\mathrm{\Phi }{m}_{\mathrm{R}}{Q}_{\mathrm{10}\mathrm{mic}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}{B}_{\mathrm{d}},\end{array}$$

where B_a→d and B_d→a denote the transition from the active to dormant microbe and from the dormant to active microbe, respectively (He et al., 2015; Wang et al., 2014). Thus, dormancy rate is affected by active and dormant biomass, soil temperature (temp), and soil moisture (θ in Φ).

The active microbial biomass (B_a) is modeled as (He et al., 2015; Wang et al., 2014)

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{B}_{\mathrm{a}}}{\mathrm{d}t}}=\text{ASSIM}\times Y_{\mathrm{g}}-m_{\mathrm{R}}{Q}_{\mathrm{10}\mathrm{mic}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}{B}_{\mathrm{a}}-B_{\mathrm{a}\to \mathrm{d}}+B_{\mathrm{d}\to \mathrm{a}}\\ & -\text{DEATH-EPROD},\end{array}\end{array}$$

where DEATH and EPROD denote microbial biomass death and enzyme production, which are modeled as proportional to active microbial biomass with constant rates r_death and r_EnzProd (Allison et al., 2010):

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathrm{DEATH}=r_{\mathrm{death}}\times B_{\mathrm{a}},\text{(11)}& {\displaystyle }\mathrm{EPROD}=r_{\mathrm{EnzProd}}\times B_{\mathrm{a}},\end{array}$$

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

The dormant microbial biomass (B_d) is modeled as (He et al., 2015; Wang et al., 2014)

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}{B}_{\mathrm{d}}}{\mathrm{d}t}}=-\mathit{\beta }{\mathrm{m}}_{\mathrm{R}}{Q}_{\mathrm{10}\mathrm{mic}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}{B}_{\mathrm{d}}+B_{\mathrm{a}\to \mathrm{d}}-B_{\mathrm{d}\to \mathrm{a}}.\end{array}$$

The soluble C pool is modeled as (He et al., 2015; Allison et al., 2010)

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}\mathrm{Soluble}\mathrm{C}}{\mathrm{d}t}}=\mathrm{DECAY}-\mathrm{ASSIM}+\mathrm{ELOSS}+\mathrm{DEATH},\end{array}$$

where DECAY represents the enzymatic decay of soil organic carbon (SOC), and ELOSS represents the loss of enzyme.

DECAY is regulated by enzyme biomass (ENZ), soil organic carbon (SOC), soil temperature, and substrate quality (He et al., 2015):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \mathrm{DECAY}=V_{\max }\times Q_{\mathrm{10}\mathrm{enz}}^{\frac{\mathrm{temp}-\mathrm{15}}{\mathrm{10}}}\times \mathrm{ENZ}\\ & \times {\displaystyle \frac{\mathrm{SOC}}{{\mathrm{Km}}_{\mathrm{uptake}}+\text{SOC}}}\times (\mathrm{120}-{\text{CN}}_{\mathrm{soil}}),\end{array}\end{array}$$

where V_max is the maximum SOC decay rate, and Km_uptake is the half-saturation constant for enzymatic decay.

ELOSS is modeled as a first-order process (Allison et al., 2010) to represent enzyme turnover:

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

where r_enzloss is the rate constant of enzyme loss.

The soil organic carbon pool (SOC) is modeled as

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{dSOC}}{\mathrm{d}t}}=\mathrm{litterfall}-\mathrm{DECAY},\end{array}$$

where litterfall is estimated as a function of vegetation carbon (Zhuang et al., 2010).

Last, the enzyme pool (ENZ) is modeled as

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

With the modification of microbial carbon dynamics by considering microbial life history traits, soil decomposition is changed since it is controlled by microbes. When microbial dormancy is considered, the number of active microbes that participate in soil decomposition is much less. The changes in soil decomposition directly influence the amount of soil respiration and further influence soil nitrogen (N) mineralization that determines soil N availability for plants, affecting gross primary production (GPP). Since both GPP and R_H can be affected by microbial dormancy, net ecosystem production (NEP) will also be affected.

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

2.3 Model parameterization and validation

The detailed description of parameters that are related to microbial dormancy can be found in He et al. (2015) (Table 1). Here we calibrated the MIC-TEM-dormancy at six representative sites with gap-filled monthly net ecosystem productivity (NEP, gC m⁻² per month) data in northern temperate and boreal regions (Table 2). Site-level climatic data and soil texture data were organized for driving the model. All site information can be found on the AmeriFlux network (Davidson et al., 2000). The results for model parameterization were presented in Fig. 2. We conducted the parameterization using a global optimization algorithm known as the SCE-UA (Shuffled complex evolution) method (Duan et al., 1994). An ensemble of 50 independent sets of parameters were performed based on prior ranges from literature (Table 1) to minimize the difference between the monthly simulated and measured NEP at the chosen sites. The cost function of the minimization is

$$\begin{array}{}\text{(18)}& \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. Except for the parameters of microbial dormancy, other parameters are derived directly from MIC-TEM (Zha and Zhuang, 2018). The optimized parameters were used for model validation and regional simulations.

Table 2 Site description and measured NEP data used to calibrate MIC-TEM-dormancy.

For model validation, we chose another six sites that contain monthly NEP data from the AmeriFlux network (Table 3). Four of these six sites were also used for parameterization (Table 2). However, we used the data of different observation periods for model validation for those overlapped sites. Moreover, we also conducted site-level validations with monthly soil respiration data from the AmeriFlux network and Fluxnet dataset. The site information was provided in Table 4. For these sites, we assumed 50 % of soil respiration was heterotrophic respiration (R_H) for forest (Hanson et al., 2000) and 60 % and 70 % of that was R_H for grassland (Wang et al., 2009) and tundra (Billings et al., 1977). Because there is a limited amount of available heterotrophic respiration (RH) data, we could not conduct a regional validation for all pixels in northern temperate and boreal regions. Instead, we extracted 61 sites providing data of average annual heterotrophic respiration from the ORNL global Soil Respiration Dataset (https://daac.ornl.gov/SOILS/guides/SRDB_V4.html, last access: 9 December 2020, Bond-Lamberty et al., 2018) for model validation. The site-level observed average annual R_H was used to compare with simulated annual R_H by MIC-TEM-dormancy and MIC-TEM. The MIC-TEM-dormancy was run at monthly time steps to keep consistent with the time step of MIC-TEM. Although microbial dynamics occur at fine temporal scales (Tang and Riley, 2014), we can still quantify the cumulative impacts of microbial dynamics on carbon and nitrogen cycling at monthly time steps by not changing the model structure.

Table 3 Site description and measured NEP data used to validate MIC-TEM-dormancy.

2.4 Spatial extrapolation

For historical simulations during the 20th century, two sets of regional simulations using MIC-TEM-dormancy and MIC-TEM at a spatial resolution of 0.5^∘ latitude × 0.5^∘ longitude were conducted. Our model simulation contains two parts: spin-up and transient simulation. A typical spin-up was conducted to get the model to a steady state for each spatial location, which will be used as initial conditions for transient simulations (McGuire et al., 1992). During spin-up procedure, cyclic forcing data were used to force the model run and repeated continuously until dynamic equilibrium was achieved at which the modeled state variables show a cyclic pattern or become constant. Specifically, this study used the monthly historical climate data from 1900 to 1940 to repeatedly drive the model for the spin-up. Before spin-up procedure, the model was initialized with default built-in carbon stocks (Raich et al., 1991). During transient simulations, the calibrated ecosystem-specific parameters were used for regional simulations. The previous dynamic equilibrium was used as an initial value for transient simulation. The historical climatic forcing data, including the monthly air temperature, precipitation, cloudiness, and atmospheric CO₂ concentrations, were organized from the Climatic Research Unit (CRU TS3.1) from the University of East Anglia (Harris et al., 2014). We also used gridded data of soil texture (Zhuang et al., 2015), elevation (Zhuang et al., 2015), and potential natural vegetation (Melillo et al., 1993) from literature. In our model, we assumed that soil texture, elevation, and potential natural vegetation data only vary spatially, not vary over time (Zhuang et al., 2015).

Figure 2 Comparison between observed and simulated NEP (gC m⁻² per month) at (a) Ivotuk (alpine tundra), (b) UCI-1964 burn site (boreal forest), (c) Howland Forest (main tower) (temperate coniferous forest), (d) Univ. of Mich. Biological Station (temperate deciduous forest), (e) KUOM turf-grass field (grassland), and (f) Atqasuk (wet tundra). Note: scales are different. Error bars represent standard errors among daily measured data in 1 month.

In addition, regional simulations over the 21st century were conducted under two Intergovernmental Panel on Climate Change (IPCC) climate scenarios (RCP2.6 and RCP8.5). The future climatic forcing data under these two climate change scenarios were derived from the HadGEM2-ES model, which is a member of CMIP5project213 (https://esgf-node.llnl.gov/search/cmip5/, last access: 9 December 2020). Then the regional estimations were obtained by summing up the gridded outputs for our study region. The positive simulated 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.

2.5 Parameter equifinality effects

Our previous studies using TEM have demonstrated that equifinality derived from site-level parameterization will affect the uncertainty in the estimation of regional carbon dynamics (Tang and Zhuang, 2008, 2009). Here equifinality refers to the fact that a number of sets of parameters result in model simulations that all match the data similarly well. To quantify this effect on our simulation uncertainty, we conducted ensemble 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).

Figure 3 Box plot of parameter posterior distributions that are obtained after ensemble inverse modeling for MIC-TEM-dormancy for all six sites: US-Ivo: Ivotuk (alpine tundra), CA-NS3: UCI-1964 burn site (boreal forest), US-Ho1: Howland Forest (temperate coniferous forest), US-UMB: Univ. of Mich. Biological Station (temperate deciduous forest), US-KUT: KUOM turf-grass field (grassland), US-Atq: Atqasuk (wet tundra).

3 Results

3.1 Inversed model parameters and model validation

Using the SCE-UA ensemble method, 50 independent sets of parameters were converged to minimize the objective function. Then the optimized parameters are calculated as the mean of these 50 sets of inversed parameters. The box plot of parameter posterior distributions reflects different ecosystem properties at these sites (Fig. 3). For instance, growth yield was higher in tundra types than in forests, meaning microorganisms in environments with higher energy limitation tend to enhance the efficiency of energy transportation. In addition, alpha, the maintenance weight, was also higher in tundra types than in forests. From the plot for the parameter beta, the ratio of dormant maintenance rate to specific maintenance rate for active biomass in tundra types is lower than that in forest types. Other microbially related parameters did not differentiate much among different vegetation types.

Figure 4 Comparison between observed and simulated NEP (gC m⁻² per month) at (a) Ivotuk (alpine tundra), (b) UCI-1964 burn site (boreal forest), (c) Howland Forest (main tower) (temperate coniferous forest), (d) Bartlett Experimental Forest (temperate deciduous forest), (e) Brookings (grassland), and (f) Atqasuk (wet tundra). Note: scales are different.

Table 4 Site description and measured R_H data used to validate the MIC-TEM-dormancy model.

After parameterization, the MIC-TEM-dormancy was validated with monthly NEP data for six representative ecosystems, and the comparisons between monthly observed NEP and simulated NEP were presented in Fig. 4. With the optimized parameters, the dormancy-based model was used to reproduce NEP to compare with the measured NEP (Table 5). The R² ranges from 0.67 for Atqasuk to 0.93 for Bartlett Experimental Forest (Table 5). Generally, our new model performs better for forest ecosystems than for tundra ecosystems. Compared with MIC-TEM, the dormancy model performs better for alpine tundra, temperate coniferous forest, and grassland. For other sites, both models show similar performance (Table 5). In addition, a set of monthly soil respiration data were selected to evaluate the estimated R_H. The comparisons between monthly observed R_H and simulated R_H from two contrasting models were conducted (Fig. 5). MIC-TEM-dormancy has higher R² and lower root-mean-square error (RMSE) (Table 6). Sixty-one sites with average annual R_H in northern temperate and boreal regions were used to further evaluate the new model performance. The dormancy model has lower intercept and slope with R² of 0.45, while R² of MIC-TEM is 0.3 (Fig. 6). These analyses indicate that the new model is more realistic in representing R_H by considering microbial dormancy.

Table 5 Model validation statistics for the dormancy model and MIC-TEM at six sites with NEP data.

Figure 5 Comparison between observed and simulated R_H (gC m⁻² per month) at (a) US-EML (alpine tundra), (b) CA-SJ2 (boreal forest), (c) US-Ho2 (temperate coniferous forest), (d) US-UMB (temperate deciduous forest), (e) US-Ro4 (grassland), and (f) RU-Che (wet tundra). Note: scales are different.

Figure 6 Linear regression between simulated and observed annual R_H (gC m⁻² yr⁻¹) for (a) MIC-TEM-dormancy and (b) MIC-TEM.

3.2 Regional carbon dynamics during the 20th century

Regional extrapolation with both models estimated a regional terrestrial ecosystem carbon sink but with different magnitudes (Fig. 7c). With optimized parameters, MIC-TEM estimated a regional carbon sink of 77.6 Pg with the interannual standard deviation of 0.21 Pg C yr⁻¹ during the 20th century. However, MIC-TEM-dormancy nearly doubles the sink at 153.5 Pg with the interannual standard deviation of 0.12 Pg C yr⁻¹ during the last century (Fig. 7c). At the end of the century, MIC-TEM estimated that NEP reaches 1.0 Pg C yr⁻¹ in comparison with MIC-TEM-dormancy estimates of 1.5 Pg C yr⁻¹ (Fig. 7c). Both models simulated similar trends for regional net primary production (NPP), R_H, and NEP (Fig. 7). Generally, they show an increasing trend in the 20th century (Fig. 7). Meanwhile, with optimized parameters, MIC-TEM-dormancy estimated NPP and R_H at 7.94 and 6.4 Pg C yr⁻¹, which are 5.8 % and 16.3 % less than the estimations from MIC-TEM, respectively (Fig. 7a and b). This pronounced difference of NEP between two models comes from the disparity between the simulated NPP and R_H with them since NEP is calculated as the difference between NPP and R_H. Without considering dormancy, MIC-TEM estimates more active microbial biomass, hence overestimating both R_H and NPP (due to higher simulated N mineralization and uptake by plants), but resulting in lower NEP than that calculated by MIC-TEM-dormancy.

Figure 7 Simulated annual net primary production (NPP, a), heterotrophic respiration (R_H, b), and net ecosystem production (NEP, c) during the 20th century by the dormancy model and MIC-TEM.

Table 6 Model validation statistics for the dormancy model and MIC-TEM at six sites with R_H data.

Temporally, both models projected higher NPP and R_H in summer than in winter (Fig. 8a and b) due to higher soil temperature and moisture (McGuire et al., 1992). Setting the R_H projection from MIC-TEM as a baseline, MIC-TEM-dormancy projected 33 % less R_H in summer (May to September) and 30 % more in winter (other months) (Fig. 8b), indicating that without dormancy, the model tends to estimate lower soil respiration due to ignorance of dormant respiration in winter but higher soil respiration due to higher active biomass in summer. NEP seasonality estimated with two models are close to each other (Fig. 8c), but the dormancy model projected slightly higher NEP in summer.

Figure 8 Regional annual seasonal pattern of simulated (a) net primary production (NPP, a), (b) heterotrophic respiration (R_H, b), and (c) net ecosystem production (NEP, c) during the 1990s from the dormancy model and MIC-TEM. The region is all land areas north of 45^∘ N.

3.3 Regional carbon dynamics during the 21st century

Under the RCP8.5 scenario, both models estimated the regional natural terrestrial ecosystems act as a carbon sink (Fig. 9). The MIC-TEM-dormancy predicted a C accumulation of 129.9 Pg by the end of this century, with the interannual standard deviation of 0.13 Pg C yr⁻¹, whereas MIC-TEM estimates a C accumulation of 79.5 Pg with the interannual standard deviation of 0.37 Pg C yr⁻¹ during the 21st century (Fig. 9). Thus, MIC-TEM-dormancy estimates an increase of 50.4 Pg regional carbon sequestration relative to MIC-TEM, with less interannual variation (Fig. 9). Under this scenario, both models predict similar temporal trends for NEP, namely increasing from the 2000s and then decreasing from the 2070s onward (Fig. 9). MIC-TEM-dormancy predicts that the carbon sink reaches 1.36 Pg C yr⁻¹ in the 2090s, which is 0.26 Pg C yr⁻¹ more than the projection of MIC-TEM. Moreover, MIC-TEM-dormancy estimated NPP and R_H at 10.2 and 8.9 Pg C yr⁻¹, which are 1.3 and 1.8 Pg C yr⁻¹ less than the estimations from MIC-TEM, respectively (Fig. 9).

Figure 9 Predicted changes in carbon fluxes: (i) NPP, (ii) R_H, and (iii) NEP for all land areas north of 45^∘ N in response to transient climate change under the RCP8.5 scenario (left panels) and RCP2.6 scenario (right panels) with the dormancy model and MIC-TEM, respectively. The decadal running mean is applied.

Under the RCP2.6 scenario, the cumulative NEP from two models diverged by 125.2 Pg C by 2100. The trajectory of interannual NEP estimated with the two models also diverged. The MIC-TEM predicted the region fluctuates between carbon sinks and sources and totally acts as a carbon source of 1.6 Pg C with the interannual standard deviation of 0.24 Pg C yr⁻¹ during the 21st century. In contrast, MIC-TEM-dormancy projected the region acts as a carbon sink of 123.6 Pg C with an interannual standard deviation of 0.1 Pg C yr⁻¹ (Fig. 9). MIC-TEM-dormancy estimates NPP and R_H at 9.9 and 8.7 Pg C yr⁻¹, which are 0.5 and 1.7 Pg C yr⁻¹ less than the estimations from MIC-TEM, respectively (Fig. 9). Moreover, simulations under the two contrasting climate scenarios (RCP2.6 and RCP8.5) exhibit a large difference of 81.1 Pg C of cumulative NEP during the 21st century by MIC-TEM, but only 6.3 Pg C of that by MIC-TEM-dormancy.

Figure 10 Regional annual seasonal pattern of simulated net primary production (NPP, a, d), heterotrophic respiration (R_H, b, e), and net ecosystem production (NEP, c, f) during the 2090s from the dormancy model and MIC-TEM under (a, b, c) the RCP2.6 scenario and (d, e, f) the RCP8.5 scenario. The region is all land areas north of 45^∘ N.

MIC-TEM-dormancy estimated higher R_H in winter, but lower R_H in summer under both future scenarios in the 2090s (Fig. 10). NPP is the same in winter with or without dormancy, and in the late summer it is higher than that without dormancy, especially in the RCP8.5 scenario. The combined flattening patterns of NPP and R_H result in different patterns for NEP. Under the RCP2.6 scenario, MIC-TEM-dormancy predicts higher NEP from June to October but lower NEP from January to April compared to MIC-TEM (Fig. 10). Under the RCP8.5 scenario, MIC-TEM-dormancy predicts higher NEP from June to September but much lower NEP in other months than MIC-TEM (Fig. 10).

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

3.4 Regional uncertainty considering equifinality effects during 20th and 21st centuries

The ensemble simulations for the 20th century are shown in Fig. 11. Given the uncertainty in parameters, MIC-TEM-dormancy predicts that the regional cumulative carbon ranges from a carbon loss of 28.2 Pg to a carbon sink of 362.1 Pg by different ensemble members, with a mean of 71.2±54.8 Pg (Fig. 11). For the 21st century, MIC-TEM-dormancy predicts that the region acts from a carbon source of 49.3 Pg C to a carbon sink of 296.5 Pg C, with a mean of 112.7±116.5 Pg under RCP2.6 scenario (Fig. 12). Under the RCP8.5 scenario, MIC-TEM-dormancy predicts that the region acts from a carbon source of 27.1 Pg C to a carbon sink of 401.3 Pg C, with a mean of 143.1±162.5 Pg (Fig. 12).

4 Discussion

Our regional simulations with two contrasting models (MIC-TEM, MIC-TEM-dormancy) indicate the regional natural terrestrial ecosystems acted as a carbon sink in past decades, which is consistent with results from other process-based models (White et al., 2000; McGuire et al., 2009; Schimel, 2013). However, the magnitudes of this sink are quite different in two models. Moreover, MIC-TEM-dormancy predicts the sink will decrease under both RCP8.5 and RCP2.6 scenarios during the 21st century, while MIC-TEM projects that the sink will increase under the RCP8.5 but change to a carbon source under the RCP2.6 scenario. Estimations based on models without dormancy could fit observations of R_H as well as estimations with dormancy, but at the cost of underestimating microbial biomass (Wang et al., 2014). Differences in predicted R_H with and without dormancy increase with temperature and with the length of the dry periods between wetting events (Salazar et al., 2018). The large difference in the two models suggests the importance of incorporating microbial dormancy effects.

Figure 12 Simulated annual net primary production (NPP, a, d), heterotrophic respiration (R_H, b, e), and net ecosystem production (NEP, c, f) under the RCP8.5 scenario (a–c) and RCP2.6 scenario (d–f) by MIC-TEM-dormancy with an ensemble of parameters. The decadal running mean is applied. The grey area represents the upper and lower bounds of simulations.

The large bias between dormancy and non-dormancy models mainly comes from two parts. First, many important microbial activities such as soil organic carbon decomposition and nutrient cycling largely depend on the active fraction of microbial communities, not total microbial biomass (Wang et al., 2014; Blagodatsky et al., 2000). However, only a small part (about 0.1 %–2 %, seldom exceed 5 %) of the total soil microbial biomass is recognized to be active under natural conditions (Blagodatsky et al., 2011; Werf and Verstraete, 1987). Thus, dormancy could be a prominent feature in soil systems (Wang et al., 2014). Without considering dormancy, the “effective” microbial biomass for soil decomposition could be overestimated, resulting in overestimation of heterotrophic respiration (He et al., 2015). He et al. (2015) predicted total soil R_H of all temperate forests (25–50^∘ N) from the dormancy model amounted to 7.28 and 8.83 Pg C yr⁻¹ from a no-dormancy model, which is 21.3 % higher than the dormancy model. Although their study region and simulation period are different from our study, the results can still be comparable. Both studies indicated that the magnitude of R_H from the no-dormancy model is higher than for dormancy models. Second, high soil respiration stimulates N mineralization in soils (Zhuang et al., 2001, 2002), making more nutrients for photosynthesis of plants (Raich et al., 1991; McGuire et al., 1995; Zhuang et al., 2015; Zha and Zhuang, 2018; Thullner et al., 2005).

Therefore, NPP will be higher due to the N enrichment from higher R_H. However, how NEP will change is still unclear. Our estimates of the northern extratropical NEP in the 1980s (1.61 Pg C yr⁻¹ with MIC-TEM-dormancy and 0.84 Pg C yr⁻¹ with MIC-TEM) are within ranges (0.6 to 2.3 Pg C yr⁻¹) reported in the literature for northern regions (Schimel et al., 2001). Moreover, our predicted time trajectory of NEP in the 21st century under the RCP2.6 scenario is very similar to the finding of White et al. (2000), indicating that NEP increases from the 2000s to the 2070s and then decreases in the 2090s. Although our dormancy model can project reasonable carbon fluxes and indicates the importance of incorporating microbial dormancy when compared with MIC-TEM (Zha and Zhuang et al., 2018), there are some other microbial traits that have not yet been considered in our model. For instance, one vital common evolutionary trait of microbes is the community shift (Wang et al., 2015) with changing environment, including warming, N fertilization, and precipitation (Treseder et al., 2011; Frey et al., 2013; Allison et al., 2009; Evans and Wallenstein, 2011). Community shift will influence microbial physiology, temperature sensitivity, and growth rates (Classen et al., 2015), which will further affect the rate of soil decomposition and other carbon dynamics (Treseder et al., 2011; Schimel and Schaeffer, 2012; Todd-Brown et al., 2011). In addition, microbial community composition was ignored in our model. We did not separate among functional microbial groups, but gather microbes into one “box”. However, microbial community composition could influence ecosystem functioning, and their variance in responses to environmental conditions could alter the prediction of the rates of decomposition of organic material (Balser et al., 2002; Fierer et al., 2007). In particular, some narrowly distributed functions can be more sensitive to microbial community composition, and these might benefit most from explicit consideration of distinguishing functional groups in ecosystem models (McGuire and Treseder, 2010; Schimel, 1995). Thus, functional dissimilarity in microbial communities can be considered the next step for model development (Strickland et al., 2009; Moorhead et al., 2006). Moreover, microbial acclimation, a mechanism of adaption to a new temperature regime, is another important trait to affect soil decomposition. Recent studies have found that the warming-induced elevated respiration of the microbial community could decrease over time because of acclimation (Melillo et al., 1993; Todd-Brown et al., 2011). This mechanism shall be factored into future soil decomposition analysis.

Except for the model limitations mentioned above, additional uncertainties may come from inadequate model parameterization and model assumptions. For example, a critical microbial parameter, carbon use efficiency (CUE), is a primary control to soil CO₂ efflux. Higher CUE indicates more microbial growth and more carbon uptake by plants, while lower CUE indicates higher soil decomposition (Manzoni et al., 2012). Theoretical and empirical studies have suggested that CUE depends on both temperature and substrate quality (Frey et al., 2013) and decreases as temperature increases and nutrient availability decreases (Manzoni et al., 2012). Our study considered the CUE sensitivity to temperature, but not nutrient availability. On the other hand, some model assumptions can also cause uncertainties. For example, we assumed that vegetation will not change during the transient simulation. However, over the past few decades in northern temperate and boreal regions, temperature increases have led to vegetation shift from one type to another (Hansen et al., 2006; White et al., 2000). The vegetation changes will affect carbon cycling in these ecosystems.

While our analysis suggests it is important to incorporate microbial dormancy dynamics into a process-based biogeochemistry model to more adequately simulate carbon dynamics in northern temperate and boreal regions, we do confront modeling dilemmas. First, our process-based models have a relatively large number of parameters, which unavoidably creates the “equifinality” problem as recognized in our previous studies for the model (e.g., Tang and Zhuang, 2008, 2009). To alleviate this problem in this analysis, we have conducted parameter ensemble simulations at both site and regional levels and presented our results with uncertainties, which could be a standard approach for process-based complex biogeochemistry modeling analyses. Second, incorporating more ecosystem processes increases the number of parameters in our model, inducing even larger uncertainties for both site-level and regional simulations. On the one hand, the more complex model to a certain degree helps capture observations; on the other hand, the model uncertainty has not been constrained or even enlarged. We highlight the need to further investigate this trade-off within the modeling research community.

5 Conclusions

This study incorporated microbial dormancy into a detailed microbe-based soil decomposition biogeochemistry model to examine the fate of large soil carbon storage in northern temperate and boreal natural terrestrial ecosystems under changing climate conditions. Regional simulations using MIC-TEM-dormancy indicated that, over the 20th century, the region is a carbon sink of 166.8±97.7 Pg. This sink could decrease to 175.9±105.4 Pg under the RCP8.5 scenario or 125.4±85.5 Pg under the RCP2.6 scenario during the 21st century. Whether considering microbial dormancy or not can cause large differences in soil decomposition estimation between two models. Meanwhile, due to available nitrogen affected by soil decomposition, net primary production is consequently influenced in these two centuries. The combined changes in soil decomposition and net primary production led to large differences in carbon budget estimation between two models. Compared with MIC-TEM, MIC-TEM-dormancy projected 75.9 Pg more C stored in the terrestrial ecosystems over the last century and 50.4 and 125.2 Pg more C under the RCP8.5 and RCP2.6 scenarios, respectively. This study highlights the importance of the representation of microbial dormancy in earth system models in order to adequately quantify the carbon dynamics of natural terrestrial ecosystems in northern temperate and boreal regions.