Modeling microbial carbon fluxes and stocks in global soils from 1901 to 2016

Soil microbes play a crucial role in the carbon (C) cycle; however, they have been overlooked in predicting the terrestrial C cycle. We applied a microbial-explicit Earth system model – the Community Land Model-Microbe (CLM-Microbe) – to investigate the dynamics of soil microbes during 1901 to 2016. The CLM-Microbe model was able to reproduce the variations of gross (GPP) and net (NPP) primary productivity, heterotrophic (HR) and soil (SR) respiration, microbial (MBC) biomass C in fungi (FBC) and bacteria (BBC) in the top 30 cm and 1 m, and dissolved (DOC) and soil organic C (SOC) in the top 30 cm and 1 m during 1901–2016. During the study period, simulated C variables increased by approximately 12 PgC yr⁻¹ for HR, 25 PgC yr⁻¹ for SR, 1.0 PgC for FBC and 0.4 PgC for BBC in 0–30 cm, and 1.2 PgC for FBC and 0.7 PgC for BBC in 0–1 m. Increases in microbial C fluxes and pools were widely found, particularly at high latitudes and in equatorial regions, but we also observed their decreases in some grids. Overall, the area-weighted averages of HR, SR, FBC, and BBC in the top 1 m were significantly correlated with those of soil moisture and soil temperature in the top 1 m. These results suggested that microbial C fluxes and pools were jointly governed by vegetation C input and soil temperature and moisture. Our simulations revealed the spatial and temporal patterns of microbial C fluxes and pools in response to environmental change, laying the foundation for an improved understanding of soil microbial roles in the global terrestrial C cycle.

1 Introduction

The atmospheric concentration of carbon dioxide (CO₂) has been drastically increased due to fossil fuel combustion and land-use change since the Industrial Revolution (IPCC, 2001, 2013; Lal, 2004, 2008). The radiative forcing caused by the CO₂ enrichment in the atmosphere has led to an increase in the global surface temperature, known as climate warming (IPCC, 2001). The increases of atmospheric CO₂ and induced warming climate have induced cascading environmental issues and impacted the carbon (C) cycle (Matson et al., 2002; Meeran et al., 2021; Soong et al., 2021).

Previous studies have assessed the effects of climate change on the global C cycle using Earth system models (ESMs) (Bonan et al., 2019; Todd-Brown et al., 2013). For example, Bonan et al. (2019) compared vegetation productivity, heterotrophic respiration, and vegetation and soil C stocks in the Community Land Model (CLM) forced by two climate reconstructions (CRUNCEPv7 and GSWP3v1) (Dirmeyer et al., 2006; Viovy, 2018). These models, however, were developed with an implicit representation of microbial processes, assuming that respired CO₂ is proportional to the soil C stock and leaving unspecified the role of microbes in decomposition processes. Given the critical role of soil microorganisms in soil biogeochemical processes and their sensitivity to environmental changes, explicit incorporation of soil microbial respiration and activities in decomposition processes into ESMs is essential to improve the prediction of global C cycling (He et al., 2021a; Wang et al., 2015, 2017; Wieder et al., 2013). Recently, researchers have applied microbial-explicit models in investigating responses of global C cycle to environmental change. For example, Wieder et al. (2015) examined the responses of soil, vegetation, and litter C pools to environmental change using the MIMICS model. Wang et al. (2017) also investigated the impacts of environmental change on enzymes, soil, and microbial biomass C pools using the TRIPLEX-MICROBE model. However, the validation of microbial biomass at coarse scales (e.g., global or biome levels) may introduce uncertainties in the model, particularly in soil microbial biomass and microbe-mediated processes, which can further affect the predicted soil C cycle in those models.

Fungi and bacteria, the two major soil microbial groups, respond differently to environmental change, and differences in their physiological traits concerning biogeochemical processes have been incorporated into the CLM-Microbe model (He et al., 2021a, b). For example, fungi decrease more than bacteria under N fertilization (Demoling et al., 2008), whereas fungi are less sensitive than bacteria to water stress (Manzoni et al., 2012). Therefore, validating fungal and bacterial biomass in the CLM-Microbe model at the grid level instead of coarse comparisons at global or biome levels may reduce uncertainties in model predictions. Changes in fungal and bacterial abundance can primarily affect terrestrial C cycling considering their distinct roles in biogeochemical processes such as the decomposition of organic materials (Bailey et al., 2002; Boer et al., 2005; Hršelová et al., 1999). Predicting changes in the spatial pattern of fungi and bacteria at the global scale and identifying their controls are essential for understanding the impacts of environmental changes on the global terrestrial C cycle.

To fill the gaps, we investigated the effects of environmental change on the global C cycle using the CLM-Microbe model. The CLM-Microbe model, mechanistically representing microbial mechanisms of soil C cycling and differentiating the physiology of two major microbial functional groups (i.e., fungi and bacteria), provides a feasible way to investigate the effects of environmental change on soil C cycling mediated by soil microbes and reveal the roles of different kingdoms of microbes on C cycling (He et al., 2021a). In this study, we aimed to investigate the effects of environmental change on the soil microbial C fluxes and pools from 1901 to 2016. We first evaluated the performance of the CLM-Microbe model in reproducing soil, vegetation, and microbial C variables, including gross (GPP) and net (NPP) primary productivity, fungal (FBC) and bacterial (BBC) biomass C in the top 30 cm and 1 m, heterotrophic (HR) and soil (SR) respiration, and dissolved organic C (DOC) and soil organic C (SOC) in the top 30 cm and 1 m. Then, we investigated the effects of environmental change on the temporal trend of microbial C fluxes and stocks including HR, SR, and FBC and BBC in the top 30 cm and top 1 m from 1901 to 2016. Finally, we investigated spatial patterns and external environmental controls of changes in those fluxes and pools from 1901 to 2016.

2 Materials and methods

2.1 Model representation of fungal and bacterial biomass

The CLM-Microbe model was built on the model framework developed by Xu et al. (2014) and the default CLM4.5 (hereafter CLM4.5) (Koven et al., 2013). It has been coupled with a microbial functional group-based methane module (Wang et al., 2019, 2022; Xu et al., 2015). Also, it has been applied to reproduce fungal and bacterial biomass dynamics and investigate fungal and bacterial macroecology and microbial seasonality impacts on soil C emission in natural ecosystems (He et al., 2021a, b, 2023). Taken together, the CLM-Microbe model has unique modules of microbe-mediated decomposition cascades and microbial functional group-mediated methane cycle, with other biogeochemical, thermal, and hydrological processes the same as the CLM4.5. The CLM-Microbe model classifies litter into three pools, i.e., litter 1 (labile), litter 2 (cellulose), and litter 3 (lignin), and soil organic matter (SOM), materials left during later stages of organic C decay, into four pools, i.e., SOM 1, SOM 2, SOM 3, and SOM 4 (low–high recalcitrance). The three litter pools and four SOM pools differ in base decomposition rate (τ), with turnover times of litter pools ranging from 20 h to 71 d and turnover times of SOM pools ranging from 14 d to 27 years (Fig. S1). Coarse woody debris (CWD) is fragmented, decomposed, and gradually transferred into litter pools and further from litter to SOM pools (Thornton et al., 2007; Koven et al., 2013). In addition to eight C pools (three litter, four SOM, and CWD pools) in the CLM4.5, we introduced dissolved organic matter (DOM) and fungal and bacterial biomass pools in the CLM-Microbe model. The code for the CLM-Microbe model has been archived at GitHub since 2015. The model version used in this study was checked out on 1 May 2021 and was archived at Xu et al. (2022). More details about the CLM-Microbe model can be found in our previous publications about the model development and model applications (Xu et al., 2014; Wang et al., 2019; He et al., 2021a, b; Wang et al., 2022; Zuo et al., 2022; He et al., 2023).

In the CLM-Microbe model, fungal and bacterial biomasses are the balance between C assimilation (C flow from the decomposition of SOM, DOM, and litter) and C loss through microbial lysis and microbial respiration. Specifically, fungi and bacteria receive C through the transitions from litter, DOM, and SOM pools; fungi and bacteria lose C through the transitions from fungal and bacterial biomass pools to DOM and SOM pools and the atmosphere. The conceptual diagram of the CLM-Microbe model and major parameters are in Fig. S1 and Table S1, respectively.

The decomposition rates of SOM, DOM, and litter are controlled by both their potential decomposition rates and environmental conditions. The decomposition processes in the CLM-Microbe model are defined following the below equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }{D}_{\mathrm{C}}=k\times r_{\mathrm{depth}}\times r_{T_{\mathrm{soil}}}\times r_{\mathrm{water}}\times r_{{\mathrm{O}}_{\mathrm{2}}}\text{(2)}& {\displaystyle }{r}_{\mathrm{depth}}=\exp (-{\displaystyle \frac{z}{z_{\mathit{\tau }}}})\text{(3)}& {\displaystyle }{r}_{T_{\mathrm{soil}}}=Q_{\mathrm{10}}^{\frac{T_{\mathrm{soil},j}-T_{\mathrm{ref}}}{\mathrm{10}}}\text{(4)}& {\displaystyle }{r}_{\mathrm{water}}=\left\{\begin{array}{cc}\mathrm{0}& for{\mathit{\phi }}_j<{\mathit{\phi }}_{\min }\\ \frac{\log \left({\mathit{\phi }}_{\min }/{\mathit{\phi }}_j\right)}{\log \left({\mathit{\phi }}_{\min }/{\mathit{\phi }}_{\max }\right)}& \text{for}{\mathit{\phi }}_{\min }\le {\mathit{\phi }}_j\le {\mathit{\phi }}_{\max }\\ \mathrm{1}& \text{for}{\mathit{\phi }}_j>{\mathit{\phi }}_{\max }\end{array}\right.\text{(5)}& {\displaystyle }\begin{array}{rl}{r}_{{\mathrm{O}}_{\mathrm{2}}}& =f_r\times \left(\mathrm{1}-f_{\mathrm{inun}}\right)\times max\left({\mathrm{O}}_{{\mathrm{2}}_{\mathrm{unsat}}},{\mathrm{O}}_{{\mathrm{2}}_{\min }}\right)\\ & +f_{\mathrm{inun}}\times max\left({\mathrm{O}}_{{\mathrm{2}}_{\mathrm{sat}}},{\mathrm{O}}_{{\mathrm{2}}_{\min }}\right),\end{array}\end{array}$$

where D_C is the rate of substrate (e.g., SOM, DOM, and litter) breakdown (in per day); k is the potential decomposition rate (in per day); $r_{{\mathrm{O}}_{\mathrm{2}}}$ represents the environmental modifier determined by soil oxygen concentration (unitless); r_depth is the environmental modifier determined by soil depth (unitless); r_water is environmental modifier determined by soil moisture (unitless); $r_{T_{\mathrm{soil}}}$ means the environmental modifier determined by soil temperature (unitless); z means soil depth (in m); z_τ is the e-folding depth for decomposition (in m); T_soil, j is soil temperature at layer j (in kelvin); T_ref is the reference temperature for decomposition (in kelvin), which is set as a kelvin temperature equals to 25 °C; Q₁₀ indicates the temperature dependence of decomposition, which is the ratio of the rate at a specific temperature to that at 10 °C lower (unitless); Ψ_j is the soil water potential in layer j (in MPa); Ψ_min is a lower limit for soil water potential control on decomposition rate (set to −10 MPa); r_water will be set as 0 if Ψ_j is lower than Ψ_min (in MPa); Ψ_max is the upper limit for soil water potential control on decomposition, which equals to the saturated soil matric potential; r_water will be set as 1 if Ψ_j is higher than Ψ_max; w_soil, j means soil water content in layer j(in MPa); f_r is the rooting fraction by soil depth (unitless); f_inun means the fraction of inundated area (unitless); O${}_{{\mathrm{2}}_{\mathrm{unsat}}}$ represents the oxygen available to that demanded by roots and aerobic microbes in unsaturated area (unitless); O${}_{{\mathrm{2}}_{\min }}$ denotes the ratio between minimum anaerobic decomposition rate and potential aerobic decomposition rate in soil (set to 0.2) (unitless); O${}_{{\mathrm{2}}_{\mathrm{sat}}}$ represents the oxygen available to that demanded by roots and aerobic microbes in saturated area (unitless); $r_{{\mathrm{O}}_{\mathrm{2}}}$ will be set as 1 in oxic conditions, while it will be estimated as the weighted average of oxygen stress in saturated and unsaturated areas in anoxic conditions (unitless).

Carbon use efficiency (CUE) of soil microbes for assimilating three litter pools in the CLM-Microbe model is determined following the equation in Sinsabaugh et al. (2013). In addition, CUE is reported to vary with temperature, showing a coefficient of −0.012 with increasing temperature (Devêvre and Horwáth, 2000). Therefore, we assumed that CUE decreased compared with the ambient thermal regime of microbes' habitats following the equation as below (Xu et al., 2014):

$$\begin{array}{}\text{(6)}& \mathrm{CUE}=\left({\mathrm{CUE}}_{\max }-{\mathrm{CUE}}_T\times \left(T-T_{\mathrm{CUEref}}\right)\right)\times {\left(M_{\mathrm{C}:\mathrm{N}}/S_{\mathrm{C}:\mathrm{N}}\right)}^{\mathrm{0.6}},\end{array}$$

where CUE is carbon use efficiency, which is defined as the growth-to-assimilation ratio for soil microbes; CUE_max is the maximum value of C use efficiency; CUE_T is the coefficient indicating the dependence of C use efficiency on temperature; T_CUEref is the reference temperature of C use efficiency, which is defined as 15 °C in the CLM-Microbe model; M_C : N means the C : N ratio of soil microbial biomass, which is defined as 8 in the CLM-Microbe model; S_C : N represents C : N ratio of the substrate (e.g., litter).

The C flow from litter and SOM pools to soil microbes will be partitioned by fungal and bacterial biomass pools based on the C : N ratio of fungal and bacterial biomass. The fraction factor quantifying bacteria C gain from litter and SOM is calculated based on the weighted average of assimilation efficiency of fungi and bacteria following the equation as below:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& fb={\displaystyle \frac{{\left(B_{\mathrm{C}:\mathrm{N}}/S_{\mathrm{C}:\mathrm{N}}\right)}^{\mathrm{0.6}}}{{\left(F_{\mathrm{C}:\mathrm{N}}/S_{\mathrm{C}:\mathrm{N}}\right)}^{\mathrm{0.6}}+{\left(F_{\mathrm{C}:\mathrm{N}}/S_{\mathrm{C}:\mathrm{N}}\right)}^{\mathrm{0.6}}}}\\ & ff=\mathrm{1}-fb,\end{array}\end{array}$$

where fb is the fraction of C flowing into bacteria; ff is the fraction of C flowing into fungi; B_C : N means the C : N ratio of BBC; F_C : N means the C : N ratio of FBC; S_C : N represents C : N ratio of substrates (e.g., litter and SOM).

Fungi and bacteria have different turnover times; hence, different lysis rate constants were adopted for fungi and bacteria in the CLM-Microbe model (He et al., 2021a). In addition, bacterial and fungal growth is highly sensitive to environmental conditions, such as soil moisture and temperature. As a result, in the CLM-Microbe model, fungal and bacterial biomass lysis is represented as the interactive effects of their lysis rate constants and environmental factors, i.e., $r_{{\mathrm{O}}_{\mathrm{2}}}$, r_water, $r_{T_{\mathrm{soil}}}$, and r_depth, as described above. Microbial respiration is widely affected by multiple abiotic and biotic factors, such as substrate concentration and availability, soil moisture, and soil temperature (Gomez-Casanovas et al., 2012; Zhang et al., 2013). Therefore, in the CLM-Microbe model, fungal and bacterial respirations are represented as the interactive effects of substrates (i.e., DOM, SOM, and litter), environmental factors (i.e., $r_{{\mathrm{O}}_{\mathrm{2}}}$, r_water, and $r_{T_{\mathrm{soil}}}$), and fraction factors quantifying C being respired by fungi and bacteria in transitions (Table S1). Fungal and bacterial biomass turnover and microbial respiration are defined following below equations:

$$\begin{array}{}\text{(8)}& {\displaystyle }L=k_M\times r_{\mathrm{depth}}\times r_{T_{\mathrm{soil}}}\times r_{\mathrm{water}}\times r_{{\mathrm{O}}_{\mathrm{2}}}\text{(9)}& {\displaystyle }R=D_{\mathrm{C}}\times f_{\mathrm{resp}},\end{array}$$

where L denotes the lysis rate of fungal and bacterial biomass (in per day); k_M is the potential turnover rate of fungal (k_fungi) or bacterial (k_bacteria) biomass (in per day); R represents the microbial respiration rate (in per day); f_resp is the fraction factor defining the proportion of C released as respiration during decomposition (unitless).

The CLM-Microbe model treats N in the same framework as CLM4.5; it fully coupled C and nitrogen (N) dynamics in land components. Net N mineralization, the inorganic N supply in the soil for plant uptake, is heavily dependent on microbial immobilization of N. Microbial immobilization of N during decomposition steps depends on C : N ratio of organic materials for decomposition, the C : N ratio of fungal and bacterial biomass, and microbial CUE. The sum of potential immobilization over all immobilization in the biogeochemistry cascade is used to estimate microbial demand of mineral N. For each time step, such microbial mineral N demand is in competition with the total plant N demand of all plant functional types (PFTs) on a soil column. Once this competition has been resolved, actual immobilization is calculated as a proportion of potential immobilization, with the same proportion applied to all immobilization steps (Thornton et al., 2007). Remaining plant N demand summed over all PFTs indicates the demand-based competition between plants and microbes for soil mineral N resource on a column. Unmet plant N demand results in C supply surplus, which is translated back to the direct downregulation of photosynthetic rate and the reduction of GPP. Unmet plant N demand can also indirectly induce the reduced allocation to new growth on light capture in plants. Such consequences imply impacts of N limitation exerted by microbial competition for mineral N on plants, which can in turn affect soil microbial community through subsequent inputs of organic matter of various qualities.

2.2 Representation of fungal- and bacterial-mediated processes by column

In the CLM-Microbe model, land surface heterogeneity was represented using a hierarchical data structure, which is adapted from CLM4.5. Each land grid cell can contain multiple land units (e.g., glacier, lake, wetland, urban, vegetated land, and cropland), and each land unit can be further divided into multiple soil/snow columns. On the vegetated land units, multiple (up to 16) PFTs distinct in physiology and structure from different climate zones (e.g., needleleaf–evergreen-tree–boreal vs. needleleaf–deciduous-tree–boreal, broadleaf–evergreen-tree–tropical vs. broadleaf–deciduous-tree–tropical, and c3–arctic-grass vs. c3-non-arctic-grass) can occupy space on the column. All vegetation fluxes and state variables were defined at the PFT level, while soil fluxes and state variables were defined at the column level.

In the CLM4.5 and early versions of the CLM-Microbe model (before January 2021), parameters related to soil processes, such as decomposition, were assumed to be homogenous across data structure levels. Our previous work suggested the differences in microbial processes among biomes (He et al., 2021a); the implicitly accounted sub-grid microbial processes may introduce uncertainties in estimating soil and microbial fluxes and state variables. Since soil flux and state variables in the CLM-Microbe model are defined at the column level, we represented the heterogeneity of microbe-mediated processes by column. Each PFT shares similar physical, phylogenetic, and phenological characteristics; we thus assigned the parameter set of microbial properties by PFT. Furthermore, we determined the microbial properties of each column by the relative weight of PFTs occupied on the column, with the parameter set of the most dominant PFT adopted to represent the microbial and soil processes (e.g., fungal and bacterial biomass turnover rate, DOM degradation rate, and fungal and bacterial C assimilation proportion from SOM, litter, and DOM) on the column.

2.3 Model forcing data

The forcing data for the CLM-Microbe model include meteorological variables such as air temperature, relative humidity, incoming solar radiation, longwave radiation, precipitation rate, surface pressure, and surface winds. In this study, we used the CRUNCEP dataset to force the CLM-Microbe model, which has been widely used to force the CLM. The CRUNCEP dataset is a combination of two existing datasets, i.e., the Climate Research Center time-series (CRU TS) dataset of 0.5° × 0.5° at a monthly scale and the National Centers for Environmental Prediction (NCEP) reanalysis dataset of 2.5° × 2.5° at 6-hourly scale. In the CRUNCEP dataset, the diurnal and daily variations of variables such as the air temperature, precipitation, humidity, solar radiation, surface pressure, downward longwave radiation, and wind speed were derived from NCEP dataset, while their monthly means are bias corrected by the CRU TS dataset. This study used the CRUNCEP dataset version 7, with a spatial resolution of 0.5° × 0.5°, spanning from 1901 to 2016, to drive the model simulation (Viovy, 2018).

In addition to the meteorological data, we forced the CLM-Microbe model using time-varying CO₂ concentration, N deposition, and aerosol concentration to estimate the C cycle change in the last century, provided by the National Center for Atmospheric Research (NCAR) for forcing the CLM offline simulations. Atmospheric N deposition during 1849–2006 with a spatial resolution of 1.25° longitude × 0.9° latitude was applied for all simulations. The CO₂ concentrations remained fixed at 1850 levels (284.7 ppm) for accelerated decomposition and final runs followed by transient historical (1849–2006) changes in the transient run. The aerosol concentration in accelerated decomposition and final runs for offline simulation was prescribed at 1850 level, while aerosol concentration with a spatial resolution of 1.25° longitude × 0.9° latitude during 1765–2005 was used in the transient simulation. The transient land use and land cover change during the historical period is based on the dataset of the UNH Transient Land Use and Land Cover Change Dataset Version 1 (LUHa.v1), covering the period of 1850–2005, which was produced by University of New Hampshire research group (Louise Chini, George Hurtt, Steve Frolking; https://luh.umd.edu/readme_LUHa_v1.shtml, last access: 10 April 2018).

2.4 Model implementation

The model implementation was carried out in three stages, with the spatial resolution of the simulations being 2.5° longitude × 1.9° latitude. First, we ran the accelerated decomposition spinup to allow the system to reach its steady state (Koven et al., 2013; Thornton and Rosenbloom, 2005). We set the model simulations to 1200 years for the accelerated decomposition phase to reach the steady state (Fig. S2). Then, we ran a final spinup of 100 years to ensure the system was ready for transient simulations during 1850–2016. For the model years of 1850–1900 in transient simulations, we cycled atmospheric forcing during 1901–1910 of the CRUNCEP dataset version 7 to force the model. Then, we used the atmospheric data during 1901 to 2016 of the CRUNCEP dataset version 7 to drive the simulation between 1901 and 2016. The CLM-Microbe model was initially parameterized for fungal and bacterial mediated processes using time-series data of fungal and bacterial biomass carbon, HR, and SR by biome. The initial setting for microbial parameters by PFT was adopted from He et al. (2021a) and He et al. (2021b). Specifically, we assigned the same microbial parameters for PFTs found in a biome as their initial setting since our previous parameterizations for microbial processes were biome-specific. For example, broadleaf–evergreen-tree–tropical and broadleaf–deciduous-tree–tropical, belonging to the biome of tropical/subtropical forests, parameter set for tropical/subtropical forests biome, were applied for PFTs of broadleaf–evergreen-tree–tropical and broadleaf–deciduous-tree–tropical in initial simulations.

2.5 Validation data

Several datasets were employed in this study for model validation. To produce realistic soil conditions in the CLM-Microbe model at the grid level, we used datasets of SOC in the top 1 m soil profile from the Harmonized World Soil Database (HWSD) at 0.05° spatial resolution archived at Oak Ridge National Laboratory (Wieder, 2014) and SOC in the top 30 cm from the Global Soil Organic C Map (GSOCmap) version 1.5 at a spatial resolution of 1 km provided by Food and Agriculture Organization of the United Nations (FAO, 2018) to validate the SOC in the top 1 m and 30 cm of the CLM-Microbe model, respectively. To guarantee the reasonability of vegetation productivity, GPP and NPP of MODIS gridded datasets with a spatial resolution of 30 s during 2000–2015 were used to compare with the simulated GPP and NPP, respectively (Zhao et al., 2005). To reproduce the soil C emission flux, SR and HR from Global Gridded 1 km Annual Soil Respiration Database (SRDB) version 3 available at Oak Ridge National Laboratory were used to validate SR and HR, respectively (Warner et al., 2019). For FBC and BBC in the top 30 cm, the dataset of FBC and BBC with a resolution of 0.5° obtained from He et al. (2020) was used to validate FBC and BBC in the top 30 cm in the CLM-Microbe model, respectively. Microbial biomass C (MBC), the sum of FBC and BBC, in the top 1 m of the CLM-Microbe model outputs was compared with Xu et al. (2013) for accuracy. The DOC in 0–30 cm and 0–1 m with a resolution of 0.5° derived from Guo et al. (2020) was compared with that in the top 30 cm and 1 m, respectively, from the CLM-Microbe output for validation. More details about the datasets used for validation can be found in Table S4. Ten-year (2000–2009) averages of simulated soil, vegetation, and microbial variables from the CLM-Microbe output were calculated to compare with those from observed datasets previously described.

To assess the efficacy of the CLM-Microbe model, the available soil and vegetation variables from the CLM4.5, including GPP, NPP, HR, SR, and SOC in the top 30 cm and 1 m, were adopted for comparison. The simulation results during 1850–2014 were forced using CRUNCEP dataset version 7, with environmental changing factors, including N deposition and rising CO_2, considered in the historical simulation. The GPP, NPP, HR, SR, and SOC in the top 30 cm and 1 m were from CLM land-only release, provided by Climate Data Gateway at the National Center for Atmospheric Research (NCAR). All variables were at a resolution of 0.9° latitude × 1.25° longitude. The temporal resolutions differed among variables, with GPP, NPP, SR, and HR being saved on a monthly scale, while SOC in 0–30 cm and 0–1 m was saved on a yearly scale. Ten-year (2000–2009) averages of the CLM4.5-simulated GPP, NPP, HR, SR, and SOC (0–30 cm and 0–1 m) were calculated to represent the long-term soil and vegetation status and for comparison with observed variables.

Since observational datasets and model simulations are of different resolutions and 0.5° is the most widely used, we used the function of linint2 in NCAR Command Language to interpolate those datasets and model outputs from their original resolutions to 0.5°. To make the maps comparable, we used the nibble and extracted by mask functions provided by ArcGIS version 10.2 (ESRI, Redlands, California, USA) to make all maps consistent in geographical boundary and missing values.

2.6 Model parameterization

Although most processes in the CLM-Microbe model were adapted from the CLM4.5, the modification of microbe-mediated decomposition cascades and microbial functional group-mediated methane cycle may reduce the applicability of default parameters in the CLM4.5. Therefore, we performed the parameterization against observational data of FBC and BBC in He et al. (2021a) and HR and SR in He et al. (2021b), with at least two sites in each biome (one for calibration and the rest for validation). Before the parameterization, we guaranteed reasonable soil and vegetation conditions by comparing the simulated NPP and SOC with observational data. To calibrate the DOC simulated in the model, we collected reported observational data from previous publications. When selecting the data, we ensured that observational DOC in the top 30 cm and 1 m was from natural biomes, and each biome included at least two sites (Tables S2–3). We calibrated the DOC in the top 30 cm and 1 m by plotting the simulated DOC in the top 30 cm and 1 m against observational data, and we found good performance of the CLM-Microbe in reproducing the observed DOC in both top 30 cm (R²=0.6, P<0.0001) and 1 m (R²=0.6, P<0.0001) (Fig. S3).

Despite the good performance of the CLM-Microbe model in the calibration phase using observational data at the site scale, minor parameter adjustments were necessary to capture variations at the global scale. We optimized the model parameters related to plant, soil, and microbial processes based on SOC in the top 30 cm from the GSOC map and in the top 1 m from the HWSD dataset, vegetation GPP and NPP from MODIS, SR and HR from SRDB, FBC and BBC (0–30 cm) in He et al. (2020), MBC in Xu et al. (2013), and DOC (0–30 cm and 0–1 m) in Guo et al. (2020). We primarily focused on parameters related to plant photosynthesis (e.g., flnr) and e-folding factor for decomposition (e.g., decomp_depth_efolding) to match the reported GPP, NPP, and SOC in the top 0–30 cm and 1 m. To calibrate the model to fit the observed FBC, BBC, and DOC, we adjusted parameters related to soil microbial (k_fungi and k_bacteria) and DOC (k_dom) turnover, microbial C assimilation efficiency (m_rf_s1m, m_ rf_ s2m, m_rf_s3m, and m_rf_s4m), the proportion of C being released as respiration (m_batm_f and m_ fatm_f), plant C allocation (froot_leaf), and the N concentration of plant tissues (leafcn and frootcn) to optimize the model simulations of FBC, BBC, MBC, DOC, SR, and HR.

2.7 Model evaluation

To evaluate the model performance in capturing the spatial variation in soil and vegetation variables, we compared GPP, NPP, HR, SR, FBC, and BBC in the top 30 cm, MBC (0–1 m) and DOC and SOC (0–30 cm and 0–1 m) reported by the observational datasets and simulated averages of these variables during 2000–2009. Due to the non-normality of those variables, Spearman's rank correlation was used to evaluate the overall model performance for those variables. The Spearman's rank correlation coefficient (r_s), measuring the strength and direction of association between two ranked variables, was calculated following the equation as below:

$$\begin{array}{}\text{(10)}& r_{\mathrm{s}}={\displaystyle \frac{\mathrm{cov}\left(R\left(x\right),R\left(y\right)\right)}{{\mathit{\sigma }}_{R\left(x\right)}{\mathit{\sigma }}_{R\left(y\right)}}},\end{array}$$

where r_s is the Spearman's rank correlation coefficient; R(x) means the rankings of variable x; R(y) indicates the rankings of variable y; cov(R(x), R(y)) is the covariance of R(x) and R(y); σ_(R(x) and σ_R(y) are the standard deviations of the rankings of variable x and y, respectively.

2.8 Statistical analysis

Due to the non-normality of simulated and observed GPP, NPP, HR, SR, FBC, and BBC in the top 30 cm, MBC (0–1 m), and DOC and SOC (0–30 cm and 0–1 m), we examined the agreement between variables simulated by the CLM-Microbe model or CLM4.5 during 2000–2009 and corresponding observed values at the grid level using Spearman's rank-order correlation. Such analyses were conducted using the function of cor.test with a method of spearman in “stats” package in R (R Core Team, 2013). The differences in soil, vegetation, litter, and microbial variables between decadal averages of 1901–1910 and 2007–2016 were examined using an independent t test, conducted with the function of t.test in “stats” package in R (R Core Team, 2013) by continents and with the function of ttest in NCAR Command Language (https://www.ncl.ucar.edu, last access: 10 April 2018) by grid.

To identify external environmental controls of soil, vegetation, litter, and microbial variables, we examined the correlations between vegetation productivity and mean annual temperature (MAT) and precipitation (MAP) and correlations of soil temperature (ST) and moisture (SM) with soil, litter, and microbial variables with respect to their area-weighted averages at the grid level from 1901 to 2016. Considering the consistent but stronger environmental influence on soil and microbial variables in the top 30 cm than in the top 1 m, only correlations between external environmental factors and soil and microbial variables in the top 1 m were assessed for whether an association exists. The correlations between external environmental factors (e.g., MAP, MAT, SM, and ST) and annual averages of GPP, NPP, HR, SR, VegC, FBC, BBC, DOC, SOC, and litter C (LitC) in the top 1 m at the global level during 1901 to 2016 were estimated using the Pearson's correlation. These statistical analyses above were performed, and relevant figures (Figs. 1–4 and 8) were plotted using “graphics” (R Core Team, 2013) and “ggcorrplot” (Kassambara and Kassambara, 2019) packages in R.

To estimate the changing rate of GPP, NPP, HR, SR, FBC, and BBC in the top 30 cm, FBC, BBC, DOC, LitC, and SOC in the top 1 m, and VegC during 1901 to 2016, we conducted linear regression models for these variables with time at the grid level, with the changing rate indicated by the slope of the regression model. In addition, correlations between environmental factors (e.g., MAT, MAP, ST, and SM) and vegetation, soil, litter, and microbial variables including GPP, NPP, HR, SR, VegC, FBC, BBC, DOC, SOC, and LitC in the top 1 m at the grid level were estimated using Pearson's correlation. Such statistical analyses were performed using NCAR Command Language (https://www.ncl.ucar.edu, last access: 10 April 2018). Relevant figures (Figs. 5–7 and 9) were produced using MATLAB version 2021b (The MathWorks, Inc.).

3 Results

3.1 Model validation and comparison with the CLM4.5

The CLM-Microbe produced comparable results with most observed data and performed comparably to or slightly better than the CLM4.5 with respect to the global C budget (Table 1), latitudinal trend (Figs. 1 and S4), and individual grid (Figs. 2 and S5). The CLM-Microbe model and the CLM4.5 overestimated the GPP by 15.7 % and 7.3 %, respectively (Table 1). However, NPP simulated by the CLM-Microbe model and the CLM4.5 was overestimated by 1.3 % and underestimated by 8.1 %, respectively. Similarly, SR was overestimated in the CLM-Microbe (15.6 %) and the CLM4.5 (4.0 %) models. HR in the CLM-Microbe and CLM4.5 models was overestimated by 1.7 % and underestimated by 4.4 %, respectively. Both the CLM-Microbe model and the CLM4.5 underestimated SOC (0–30 cm) by 8.5 % and 22.4 %, respectively, while SOC (0–1 m) in the CLM-Microbe model and the CLM4.5 was overestimated by 32.4 % and underestimated by 21.4 %, respectively. The FBC, BBC, MBC, and DOC, only available in the CLM-Microbe model, were better predicted in the top 30 cm than 1 m. The simulated FBC, BBC, and DOC in the top 30 cm were underestimated by 3.3 % and overestimated by 26.7 % and 24.9 %, respectively, while MBC and DOC in the top 1 m were overestimated by 69.5 % and 75.0 %, respectively.

Table 1 Annual flux of GPP, NPP, HR, and SR and carbon stocks of FBC in the top 30 cm, BBC in the top 30 cm, MBC (0–1 m), DOC in the top 30 cm, SOC in the top 30 cm, and SOC in the top 1 m by observed datasets and by simulations of the CLM-Microbe model and CLM4.5 at the global scale.

GPP: gross primary productivity; NPP: net primary productivity; HR: heterotrophic respiration; SR: soil respiration; DOC: dissolved organic carbon; SOC: soil organic carbon; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; MBC: microbial biomass carbon. NA: not available. The SOC (0–1 m) data are from the Harmonized World Soil Database (HWSD, https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1247, last access: 10 April 2018); the SOC (0–30 cm) data are from the Global Soil Organic Carbon Map (GSOCmap) version 1.5; GPP and NPP data are from MODIS gridded datasets (http://files.ntsg.umt.edu/data/NTSG_Products/, last access: 10 April 2018); the SR and HR data are from Global Gridded 1 km Annual Soil Respiration Database (SRDB) version 3 (https://daac.ornl.gov/CMS/guides/CMS_Global_Soil_Respiration.html, last access: 10 April 2018); the FBC and BBC in the top 30 cm are from He et al. (2020); MBC (0–1 m) is compared with Xu et al. (2013); the DOC (0–30 cm and 0–1 m) is derived from Guo et al. (2020). Output of the CLM-Microbe model during 2000–2009 (decadal average) is used to compare with observational data.

The CLM-Microbe model can reasonably reproduce the latitudinal trends of vegetation, soil, and microbial variables, with the model performance varied among variables and along soil depth (Fig. 1). The latitudinal trends of both GPP and NPP in the CLM-Microbe model agreed with observed data with a slight overestimation at northern latitudes and in equatorial regions, but NPP was slightly underestimated in the Southern Hemisphere (Fig. 1a and b). Both HR and SR simulated by the CLM-Microbe model agreed well with observed data in the Southern Hemisphere but were overestimated in equatorial regions and at mid-high latitudes and underestimated at low latitudes in the Northern Hemisphere (Fig. 1c and d). Similar latitudinal trends of HR and SR were also observed in the CLM4.5 simulation (Fig. S4c and d). Soil C pools showed similar latitudinal patterns across soil depths (Fig. 1e–k). Specifically, DOC (0–30 cm and 0–1 m) was overestimated in equatorial regions but underestimated in northern temperate regions (Fig. 1e–f). Meanwhile, the CLM-Microbe model overestimated SOC (0–30 cm and 0–1 m) in equatorial and northern high-latitude regions but underestimated in northern mid-latitude regions (Fig. 1g and h). As opposed to the CLM-Microbe model, the CLM4.5 consistently underestimated SOC (0–30 cm and 0–1 m) along latitudes, except for SOC (0–1 m) at latitudes of > 60° N. Similarly, both FBC and BBC in the top 30 cm were overestimated in equatorial regions and at northern high latitudes but underestimated in northern mid-latitude regions (Fig. 1i and j). Overall, FBC (0–30 cm) at southern latitudes was well predicted by the CLM-Microbe model, but BBC (0–30 cm) in that region was underestimated, while MBC (0–1 m) was overestimated across latitudinal gradients (Fig. 1k).

Figure 1 Latitudinal comparison between observed (black line) and the CLM-Microbe-simulated (red line) (a) GPP, (b) NPP, (c) HR, (d) SR, and (e) DOC in the top 30 cm, (f) DOC in the top 1 m, (g) SOC in the top 30 cm, (h) SOC in the top 1 m, (i) FBC in the top 30 cm, (j) BBC in the top 30 cm, and (k) MBC (0–1 m). GPP: gross primary productivity; NPP: net primary productivity; HR: heterotrophic respiration; SR: soil respiration; DOC: dissolved organic carbon; SOC: soil organic carbon; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; MBC: microbial biomass carbon.

At the grid scale, the simulated values of vegetation, soil, and microbial variables with the CLM-Microbe model were significantly consistent with the observational results (P<0.05; Fig. 2). The CLM4.5 also indicated significant consistency between simulated and observed vegetation and soil variables (P<0.05; Fig. 2). Overall, the CLM-Microbe model and CLM4.5 performed well at simulating GPP and NPP compared to simulating soil and/or microbial variables. The simulated GPP (r_s=0.91) and NPP (r_s=0.86) with the CLM-Microbe model were significantly and positively related to their observed values (Fig. 2a and b). The GPP (r_s=0.88) and NPP (r_s=0.82) in the CLM4.5 were also significantly and positively associated with observed values (Fig. S5a and b). The SR tended to be better predicted than HR in both the CLM-Microbe model (r_s=0.70 for SR vs. r_s=0.68 for HR) and the CLM4.5 (r_s=0.68 for SR vs. r_s=0.64 for HR) (Fig. 2c and d and Fig. S5c and d). The DOC in 0–1 m (r_s=0.36) was slightly better reproduced than in 0–30 cm (r_s=0.34) in the CLM-Microbe model (Fig. 2e and f), while both the CLM-Microbe model (r_s=0.68 for 0–30 cm vs. r_s=0.63 for 0–1 m) and CLM4.5 (r_s=0.63 for 0–30 cm vs. r_s=0.59 for 0–1 m) performed better at simulating SOC in the top 30 cm than in the top 1 m (Figs. 2g and h, S5e and f). Similarly, the CLM-Microbe model performed better in simulating FBC and BBC in the top 30 cm than MBC in the top 1 m (r_s=0.43) (Fig. 2i–k). In addition, BBC (r_s=0.53) was better reproduced than FBC (r_s=0.46) in the top 30 cm.

Figure 2 Grid-by-grid comparison between observed and the CLM-Microbe-simulated (a) GPP, (b) NPP, (c) HR, (d) SR, and (e) DOC in the top 30 cm, (f) DOC in the top 1 m, (g) SOC in the top 30 cm, (h) SOC in the top 1 m, (i) FBC in the top 30 cm, (j) BBC in the top 30 cm, and (k) MBC (0–1 m). Red lines are 1:1 line. GPP: gross primary productivity; NPP: net primary productivity; HR: heterotrophic respiration; SR: soil respiration; DOC: dissolved organic carbon; SOC: soil organic carbon; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; MBC: microbial biomass carbon.

3.2 Carbon fluxes and pools associated with soil microbes

The HR and SR displayed increasing trends from 1901 to 2016 (Fig. 3a–b), with different magnitudes among variables. By 2016, the increase of SR (25 PgC yr⁻¹) was about twice that of HR (12 PgC yr⁻¹). Their increasing rates showed variations with time. We observed a relatively modest increase in HR and SR during 1901–1980, whereas their increases were more rapid from 1981–2016. In addition, microbial C pools increased from 1901 to 2016 despite the year-to-year variability (Fig. 3c–f). The FBC and BBC in the top 30 cm and FBC and BBC in the top 1 m increased by about 1.0, 0.4, 1.2, and 0.7 PgC, respectively, from 1901 to 2016. However, the temporal trends of those variables varied during 1901 to 2016. The FBC (0-30 cm and 0–1 m) decreased from 1901–1940 and increased after 1940 (Fig. 3c–d). The BBC (0–30 cm and 0–1 m) exhibited little change during 1901–1940 but increased rapidly during 1941–2016 (Fig. 3e–f).

Figure 3 Evolution of annual carbon flux of area-weighted (a) HR and (b) SR and carbon stock of (c) FBC in the top 30 cm, (d) FBC in the top 1 m, (e) BBC in the top 30 cm, and (f) BBC in the top 1 m simulated by the CLM-Microbe model since 1901. The baseline was the 10-year average of corresponding variables during 1901–1910. HR: heterotrophic respiration; SR: soil respiration; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; MBC: microbial biomass carbon.

3.3 Spatial pattern of microbial carbon fluxes and stocks

Compared with 1901–1910, HR and SR increased across latitudinal gradients in 2007–2016 (Fig. 4a–b). However, the magnitude of the increase differed among latitudinal gradients. Specifically, increases in HR and SR were larger and more prominent at northern latitudes and equatorial regions than at southern latitudes. Similar to C fluxes, microbial C pools increased across latitudinal gradients from 1901–1910 to 2007–2016 (Fig. 4c–f). Overall, FBC and BBC in the top 30 cm and top 1 m showed a small but different extent of increase across latitudinal gradients. Specifically, the increases were more prominent at northern high latitudes and equatorial regions than at other latitudes.

Figure 4 Latitudinal gradients of the CLM-Microbe model simulated 10-year averages of (a) HR, (b) SR, and (c) FBC in the top 30 cm, (d) FBC in the top 1 m, (e) BBC in the top 30 cm, and (f) BBC in the top 1 m during 1901–1910 and 2007–2016. HR: heterotrophic respiration; SR: soil respiration; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; MBC: microbial biomass carbon.

Across the globe, HR and SR showed similar spatial patterns, and increases in most grids across the globe were statistically significant (P<0.05; Fig. 5a, b, d, and e). Correspondingly, we observed positive relative change in most areas from 1901–1910 to 2007–2016 (Fig. 5c and f). However, we also observed decreases in HR and SR in the grids of South Asia. The HR and SR displayed similar spatial patterns of changing rates (Fig. 9). At the global scale, FBC and BBC in the top 30 cm and 1 m showed similar spatial patterns and widely increased from 1901–1910 to 2007–2016 (Fig. 6a and b, d and e, g and h, and j and k). Correspondingly, we observed positive relative changes in FBC and BBC in the top 30 cm and 1 m in those regions from 1901–1910 to 2007–2016 (Fig. 6c, f, i, and l).

Figure 5 Spatial distributions of decadal averages of (a–b) HR and (d–e) SR during (a, d) 1901–1910 and (b, e) 2007–2016 and relative changes in (c) HR and (f) SR by 2007–2016 relative to 1901–1910. HR: heterotrophic respiration; SR: soil respiration. Black dot in each grid indicates significant changes (P<0.05).

Consistent with the spatial patterns of absolute and relative changes, we widely observed significant and positive changing rates of HR and SR from 1901 to 2016 (Fig. 7a–b; P<0.05). However, we also found significant negative changing rates of HR and SR in grids of South Asia (P<0.05). Similarly, increasing temporal trends of FBC and BBC in the top 30 cm and 1 m were widely observed across the globe (Fig. 7c–f). However, we also observed decreases of those variables in South Asia (Figs. 5–7). In addition, we observed decreases of FBC and BBC in the top 30 cm and 1 m in grids of central North America (Figs. 6 and 7c–f).

Figure 6 Spatial distributions of decadal averages of (a–b) FBC in the top 30 cm, (d–e) FBC in the top 1 m, (g–h) BBC in the top 30 cm, and (j–k) BBC in the top 1 m during (a, d, g, j) 1901–1910 and (b, e, h, k) 2007–2016 and relative changes in (c) FBC in the top 30 cm, (f) FBC in the top 1 m, (i) BBC in the top 30 cm, and (l) BBC in the top 1 m by 2007–2016 relative to 1901–1910. FBC: fungal biomass carbon; BBC: bacterial biomass carbon. Black dot in each grid indicates significant changes (P<0.05).

Figure 7 Changing rates of the CLM-Microbe model simulated (a) HR, (b) SR, and (c) FBC in the top 30 cm, (d) FBC in the top 1 m, (e) BBC in the top 30 cm, and (f) BBC in the top 1 m from 1901 to 2016. HR: heterotrophic respiration; SR: soil respiration; FBC: fungal biomass carbon; BBC: bacterial biomass carbon. Black dot in each grid indicates significant regression (P<0.05).

3.4 External environmental controls on C cycling

The area-weighted average of HR, SR, and FBC and BBC in the top 1 m were significantly correlated with that of ST and SM in the top 1 m (P<0.05; Fig. 8). However, the strengths of correlations depended on both environmental controls (ST and SM in the top 1 m) and variables (HR, SR, and FBC and BBC in the top 1 m). For example, correlations of HR and SR with ST and SM in the top 1 m were of the same magnitude, while the FBC and BBC were more strongly correlated with ST than with SM in the top 1 m.

Across the globe, microbial carbon fluxes and stocks were more widely and positively correlated with ST than with SM in the top 1 m (Fig. 9). Correlations of ST (0–1 m) with HR and SR were similar in spatial patterns. We widely observed significant positive correlations of HR and SR with ST (0–1 m) (P<0.05; Fig. 9a and c). On the other hand, negative correlations of HR and SR with ST (0–1 m) were found in South Asia, southeast North America, central North America, central Africa, and central and northern Australia/Oceania. Correlations of FBC and BBC in the top 1 m displayed similar spatial patterns. We found significant and positive correlations of FBC and BBC with ST in the top 1 m in most grids across the globe (P<0.05; Fig. 9b, d, f, and h). However, we also found some grids with negative correlations in central North America, Europe, Asia, South America, Africa, and Australia/Oceania. Correlations of HR, SR, and FBC and BBC in the top 1 m were similar in spatial patterns, with significant and positive correlations widely observed (P<0.05). But we also observed negative correlations at middle and low latitudes in North America, Europe, and Asia, the east coast of South America and Africa, and southern Australia/Oceania.

Figure 8 Heatmap showing Pearson’s correlation between the CLM-Microbe model simulated microbial carbon fluxes and pools (HR, SR, FBC in the top 1 m, and BBC in the top 1 m) and soil environments (SM and ST in the top 1 m) from 1901 to 2016. GPP: gross primary productivity; NPP: net primary productivity; HR: heterotrophic respiration; SR: soil respiration; DOC: dissolved organic carbon; SOC: soil organic carbon; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; VegC: vegetation carbon; LitC: litter carbon; MAT: mean annual temperature; MAP: mean annual precipitation; ST: soil temperature; SM: soil moisture. Black asterisks indicate significant correlations (P<0.05).

Figure 9 Pearson's correlation between the CLM-Microbe model simulated (a–b) HR, (c–d) SR, and (e–f) FBC in the top 1 m, (g–h) BBC in the top 1 m, and (a, c, e, g) ST and (b, d, f, h) SM in the top 1 m from 1901 to 2016. Black dot in each grid indicates significant correlation (P<0.05). HR: heterotrophic respiration; SR: soil respiration; FBC: fungal biomass carbon; BBC: bacterial biomass carbon; ST: soil temperature; SM: soil moisture.

4 Discussion

4.1 Comparison with previous studies

The latitudinal trends and grid-level distribution of GPP, NPP, HR, SR, FBC, and BBC in the top 30 cm and FBC, BBC, DOC, and SOC in the top 1 m were well-reproduced in the CLM-Microbe model (Figs. 1 and 2). The CLM-Microbe model performed better than or comparable to the CLM4.5 in simulating the spatial distribution of vegetation, soil, and microbial variables (Figs. 2 and S5). In line with our results, multiple models captured the spatial variation of GPP, NPP, HR, and SR (Delire et al., 2020; Kim et al., 2019; Wiltshire et al., 2021; Zheng et al., 2020). Wieder et al. (2015) reported a high spatial correlation (r=0.46) of SOC (0–1 m) between MIMICS outputs and HWSD. In addition, Wang et al. (2017) observed the high consistency in SOC (0–1 m) (R²=0.96; P<0.01) between the TRIPLEX-MICROBE model and HWSD by vegetation type. Huang et al. (2021) also found good performance of the ORCHIMIC v2.0 in reproducing SOC by comparing the simulated values with multiple SOC datasets. The well-developed plant physiology and environmental controls in the model may explain their good performance in simulating vegetation and soil processes (Flato, 2011; Mathieu and O'Neill, 2008). However, the latitudinal trends and grid-level distribution of DOC, SOC, and MBC (sum of FBC and BBC) in the top 1 m were relatively worse reproduced than those in the top 30 cm (panels e–k of Figs. 1 and 2, panels e–f of Figs. S4 and S5), indicating that the vertical distribution of processes related to decomposition, microbial turnover, and plant C input needs further improvements. Although parameters classifying the active decomposition depth and biological function to perturbation were defined in the CLM-Microbe model, the gradual change of microbial turnover and activity defined along the soil profile may need to be improved in future models (Preusser et al., 2019; Zhu et al., 2021). In addition, processes or parameters related to the active layer for decomposition and perturbation caused by biological (e.g., nematode) and abiotic (e.g., drying and rewetting) activities can cause uncertainties in the vertical distribution of C cycle, which needs further efforts and attention in model development (Ettema and Wardle, 2002; Gabet et al., 2003; Kuzyakov and Blagodatskaya, 2015; Schimel, 2018).

We estimated global annual averages of 129.5, 56.5, 99.8, and 49.8 PgC yr⁻¹ for GPP, NPP, HR, and SR, respectively (Table 1). Consistent with our results, previous studies reported similar values of GPP, NPP, HR, and SR (Cramer et al., 1999; Hashimoto et al., 2015; Huang et al., 2020; Lu et al., 2021; Nemani et al., 2003; Zhao et al., 2017; Zheng et al., 2020). The consistent simulations and reasonable estimations of GPP, NPP, HR, and SR across models may indicate the convergent plant physiology among models and well-defined soil and microbial processes in the CLM-Microbe model. In addition, compared with observed data, the CLM-Microbe model produced more consistent NPP and HR but overestimated GPP and SR (Table 1). The overestimation of GPP and SR may be due to the lower ecosystem-scale CUE in the CLM-Microbe model. The vegetation physiology module in the CLM-Microbe model is adapted from CLM4.5. The ecosystem-scale CUEs between the CLM-Microbe model (0.44) and CLM4.5 (0.43) were comparable but lower than in MODIS (0.5). Correspondingly, we observed a higher contribution of roots to total SR in the CLM-Microbe model (0.5) and the CLM4.5 (0.48) than in the observed SRDB dataset (0.43). Therefore, the well-simulated NPP and HR but higher predicted GPP and SR in the CLM-Microbe model were attributed to the low ecosystem-scale CUE. Increasing ecosystem-scale CUE in the CLM-Microbe model will improve the modeling performance of GPP and SR in model development.

The CLM-Microbe model can reasonably predict FBC, BBC, and DOC in the top 30 cm well globally, indicating the well-represented microbial processes in surface soils (Table 1). However, MBC and DOC in the top 1 m were vastly overestimated, with MBC and DOC in the top 1 m overestimated by 69.5 % and 75.0 %, respectively. Inconsistent with our results, previous studies suggested the underestimation of MBC (0–1 m) in their models. For example, Wang et al. (2017) estimated the global MBC as 21 PgC in the TRIPLEX-MICROBE model. Wieder et al. (2015) suggested the steady-state MBC (0–1 m) of 16.3 Pg in the MIMICS. The relatively poor performance of the CLM-Microbe model in simulating DOC and MBC in the top 1 m and the discrepancy in simulated MBC (0–1 m) among studies may result from three aspects. First, the hydrologically active layer in the CLM-Microbe model may not be sufficient to define soil microbial processes along soil profiles. We observed better performance of the CLM-Microbe model in simulating FBC, BBC, and DOC in the top 30 cm relative to MBC and DOC in the top 1 m, indicating that the representation of microbial and soil processes along soil profiles may need improvements. Second, the difference in calibration for MBC may cause a discrepancy between studies. The SOC in Wieder et al. (2015) was calibrated to observed data but not MBC. Wang et al. (2017) calibrated the MBC (0–1 m) in the TRIPLEX-MICROBE by vegetation types, while we calibrated both MBC and SOC in 0–30 cm and 0–1 m by grid in the CLM-Microbe model. The differences in variables and depths calibrated between studies can partly explain the discrepancy. Third, the difference in simulated vegetation, litter, and soil C pools among studies can result in the discrepancy. Vegetation C as litter and volatile organic compounds, DOC, and SOC are the C source for microbial C assimilation through decomposition (Fig. S1). Consequently, the overestimation of SOC and DOC can partly explain the overestimation of MBC in the top 1 m (Table 1).

The CLM-Microbe model indicated an underestimation of 8.5 % for SOC (0–30 cm) and an overestimation of 32 % for SOC (0–1 m) when comparing with observed data (Table 1). Compared with the CLM4.5, the CLM-Microbe predicted larger stocks of SOC (0–30 cm and 0–1 m). Previous studies suggest large variations in simulated SOC (0–1 m) among models. For example, Todd-Brown et al. (2013) reported the SOC (0–1 m) stock ranging from 510 to 3040 PgC among 11 CMIP5 ESMs. The TRIPLEX-MICROBE modeled the global SOC (0–1 m) stock as 1195 PgC (Wang et al., 2017). Wieder et al. (2015) documented the steady-state SOC pool in the MIMICS as 1530 PgC. Delire et al. (2020) reported the SOC (0-1 m) as 1611 and 1520 PgC in the new (ISBA_bgc6) and old (ISBA_bgc5) versions, respectively, of ISBA-CTRIP. Given the wide range (510 to 3040 PgC) of simulated SOC (0–1 m) in models, the CLM-Microbe model thus predicted reasonable SOC stocks.

4.2 Temporal trends of carbon fluxes and stocks of soil microbes

The area-weighted average of HR and SR in the CLM-Microbe model increased by 12 and 25 PgC yr⁻¹, respectively, from 1901 to 2016 (Fig. 3). Consistent with our findings, Bonan et al. (2019) observed the increase of about 8 PgC yr⁻¹ of HR from 1850 to 2014 in the CLM4.5. The global increasing rate of SR was estimated as 0.04–0.14 PgC yr⁻¹ by Huang et al. (2020). The rising ST (0–1 m) may explain the observed increase in HR considering the positive relationship between ST (0–1 m) and HR (Figs. 8 and S6c). The increase in HR can partly explain the rising SR from 1901 to 2016 given its critical contribution to SR. In addition to HR, the increase in root respiration due to increasing C availability and rising temperature accounted for a crucial proportion of the SR increase (Bond-Lamberty and Thomson, 2010; Hashimoto et al., 2015; Piñeiro et al., 2017; Zhou et al., 2016). We observed increases of GPP and NPP associated with environmental changes such as increasing N deposition and rising CO₂ concentration and temperature (Dusenge et al., 2019; Piñeiro et al., 2017). Evidenced by increasing VegC and LitC, indicating the C stock of vegetation biomass and C loss of vegetation biomass, respectively, C input from plants increased during the historical period at the global scale (Fig. S8e–f). Therefore, increases in vegetation C sequestration, together with rising ST and SM in the top 1 m, enhanced HR and SR in the last century.

The area-weighted FBC and BBC in 0–30 cm increased by 1.0 and 0.4 PgC and those in 0–1 m increased by 1.2 and 0.7 PgC, respectively, from 1901 to 2016 in the CLM-Microbe model (Fig. 3c–f). Soil microbes are sensitive to environmental change, and rising temperature was reported to induce lower microbial biomass due to the negative impacts of temperature on microbial biomass maintenance through facilitating microbial turnover (Joergensen et al., 1990; He and Xu, 2021). We observed an increasing trend of ST (0–1 m) from 1901 to 2016 (Fig. S6c), indicating negative impacts of temperature on FBC and BBC in 0–30 cm and 0–1 m. In addition to temperature, microbial biomass is influenced by substrate and water availability, with significant and positive effects of SOC and MAP recorded on microbial biomass (Chen et al., 2022). Litter, SOM, and DOC are three C sources for soil microbes in the CLM-Microbe model (“Materials and methods”; Fig. S1). We observed increases of DOC (2.4 PgC), LitC (4 PgC), and SOC (34 PgC) in the top 1 m from 1901 to 2016, indicating more C available for soil fungi and bacteria during the historical period (Fig. S8c, d, and f). In addition, MAP and SM (0–1 m) increased from 1901 to 2016 (Fig. S6b and d). Therefore, the increasing substrates (DOC, LitC, and SOC) and soil water availability can explain the increase in FBC and BBC in the CLM-Microbe model.

The annual averages of microbial C fluxes (HR and SR) and pools (FBC and BBC in the top 30 cm and 1 m) have exhibited more rapid increases since 1980 (Fig. 3). Concurrently, we observed a more rapid increase in MAT, MAP, and ST and SM in the top 1 m since 1980 (Fig. S6). In line with this study, Cheng et al. (2017) analyzed SM simulations during historical (1920–2005) and future (2006–2080) periods in the CESM from CIMP5; they also found 1980 as a transition for a subsequent increase of variation during 1920–2005, indicating more rapid changes in SM after 1980. We observed significant correlations of HR, SR, FBC (0–1 m), and BBC (0–1 m) with ST and SM in the top 1 m (Fig. 8). Therefore, more rapid increases in MAT, MAP, and ST and SM in the top 1 m after 1980 may explain the more rapid increases of such variables since 1980.

4.3 Changes in microbial carbon fluxes and stocks over the space and their controls

The HR and SR showed an increase across latitudinal gradients and the globe in the study period (Figs. 4a and b and 5). Consistent with our findings, Huang et al. (2020) observed a globally significant increase in SR, particularly in boreal and tropical regions (e.g., northern Asia, central South America, and central and southern Africa), from 2000–2014. Bond-Lamberty et al. (2018) also observed an increase in HR in multiple biomes during 2000–2015. In addition, we observed similar spatial patterns of increases (e.g., higher increases at northern high latitudes and in equatorial regions) in HR and SR with those of GPP and NPP (Fig. S10). These results indicated that soil C fluxes largely depended on vegetation productivity, which can enhance soil C fluxes due to high C allocation to belowground (Pendall et al., 2004; Prescott et al., 2020). In addition, soil C fluxes can be further increased due to facilitated decomposition in a warming world (Noh et al., 2017; Zhou et al., 2007). Temperature and water availability have a profound influence on root respiration and HR (Bond-Lamberty and Thomson, 2010; Hashimoto et al., 2015; Sinsabaugh et al., 2016). We also found significant correlations of HR and SR with ST and SM in the top 1 m (Fig. 9). The increases in HR and SR can be explained by the increases in SM and ST in the top 1 m, considering their significant correlations (Fig. S6c and d). However, we also observed decreases in HR and SR in South Asia (Fig. 5). Vegetation C fixation is the major C source for ecosystems; the decreases of GPP and NPP in such regions can largely explain the decrease in HR and SR (Fig. S10).

The FBC and BBC in the top 30 cm and 1 m increased across latitudes during 2007–2016 compared with 1901–1910 (Fig. S4c–f). In addition, the FBC and BBC in the top 30 cm and 1 m widely increased across the globe (Figs. 6 and 7c–f). Vegetation is the major C source for soil microbes in terrestrial ecosystems, determining the total amount of C available for microbes by regulating microbial C assimilation through SOC, DOC, and litter (Schimel, 1995; Vance and Chapin, 2001; Xu et al., 2014). The spatial patterns of GPP and NPP change could explain the wide increases of FBC and BBC in the top 30 cm and 1 m in such areas, as well as their larger increases at high latitudes and in equatorial regions (Figs. S10 and S12a–b). However, we also found slight decreases of FBC and BBC in the top 1 m in regions such as southern Australia/Oceania (Figs. 6 and S7c–f). As a critical C source for soil microbes in the CLM-Microbe model, the decrease in DOC (0–1 m) may explain the widespread decrease in FBC and BBC of 30 cm and 1 m in southern Australia/Oceania (Figs. S1, S11a–c, and S12c). Meanwhile, in the top 1 m, FBC and BBC were negatively correlated with ST and positively correlated to SM in southern Australia/Oceania (Fig. 9e–h). The increase in ST and SM in the top 1 m can explain the decrease in FBC and BBC in 0–30 cm and 0–1 m in southern Australia/Oceania (Fig. S7c and d). In addition, we found decreases of FBC and BBC in the top 30 cm and 1 m in South Asia and central North America (Figs. 6 and 7c–f). Since vegetation productivity is the primary C source for terrestrial ecosystems, decreases in vegetation C input as GPP and NPP can explain their decreases in South Asia (Figs. S10 and S12a and b). In addition, microbial activities are affected by soil temperature and water availability. We observed negative correlations of FBC and BBC in the top 1 m with SM (0–1 m) in South Asia and with ST (0–1 m) in central North America. Increases of SM (0–1 m) and ST (0–1 m) in corresponding regions from 1901–1910 to 2007–2016 may contribute to decreases in FBC and BBC in the top 30 cm and 1 m in such areas (Figs. 9e–h and S7c and d).

Correlations of microbial C fluxes (HR and SR) with SM and ST in the top 1 m varied across space (Fig. 9a–d). Specifically, the association between respiration fluxes, HR and SR, and soil climatic factors (ST and SM) in the top 1 m were positive in the majority of land area (P<0.05; Fig. 9a–d) but were negative in central North America and South America, South Asia, central Africa, and central and northern Australia/Oceania. Studies found positive effects of rising temperature and increasing water availability on microbial activities (Nyberg and Hovenden, 2020; Tecon and Or, 2017); therefore, the widely increases of ST and SM in the top 1 m can explain the increases of HR at the global scale (Fig. S7c and d). In addition to soil moisture and temperature, HR is positively related to substrates and microbial biomass (Wei et al., 2015). Therefore, the reduction in HR in South Asia can be explained by decreases in DOC, LitC, SOC, FBC, and BBC in the top 1 m, while decreasing SM (0–1 m) and ST (0–1 m) may contribute to the reduced HR in central Africa and central and northern Australia/Oceania, respectively (Figs. 6, 7c–f, S11–12, and S7c–d). Since HR contributed over 50 % to SR, factors determining correlations of HR with SM and ST in the top 1 m were expected to be responsible for those of SR across space. In addition, the reduction in root respiration would enhance the negative correlations of HR with SM and ST in top 1 m in South Asia considering its decreasing vegetation productivity (Figs. S10 and S12a–b). Water availability impacts on microbial activities are affected by other factors such as temperature and substrate availability (Moyano et al., 2013; Tecon and Or, 2017). We observed decreasing SM (0–1 m) at middle and low latitudes in North America, Europe, and Asia, the east coast of South America and Africa, and southern Australia/Oceania (Fig. S7d). However, HR and SR and FBC and BBC in the top 1 m widely increased due to the positive effects of temperature and substrate availability (Figs. 5–7). Therefore, negative correlations of HR, SR, and FBC and BBC with SM in the top 1 m at middle and low latitudes in North America, Europe, Asia, the east coast of South America and Africa, and southern Australia/Oceania resulted from the facilitating effects of temperature and substrate availability on HR, SR, and FBC and BBC in the top 1 m even with decreasing SM (0–1 m).

Microbial C pools showed varied correlations with SM and ST in the top 1 m across space (Fig. 9). Specifically, FBC and BBC mostly showed significant positive correlations with ST and SM in the top 1 m (P<0.05; Fig. 9e–h). But we also observed negative correlations of FBC and BBC in the top 1 m with ST (0–1 m) in central North America and South America and with SM (0–1 m) at middle and low latitudes in North America, Europe, Asia, the east coast of South America and Africa, and southern Australia/Oceania. Temperature negatively impacts microbial biomass due to enhanced microbial turnover with rising temperature, but such a relationship depends on environmental conditions (Joergensen et al., 1990; He and Xu, 2021; Yuste et al., 2007). Water availability enhances microbial biomass. The positive correlations of FBC and BBC with SM in the top 1 m at middle and high latitudes in Asia, Europe, and North America contributed to positive correlations of FBC and BBC with ST in the top 1 m in such regions. In addition, substrate availability is critical for soil microbes. Increases in DOC, SOC, and LitC in the top 1 m can explain the widely found positive correlations of FBC and BBC with ST (0–1 m) across space (Figs. S11–12). Water can affect soil microbes in distinct ways, affected by other factors such as temperature and substrate availability (Moyano et al., 2013; Tecon and Or, 2017). We observed decreasing SM (0–1 m) at middle and low latitudes in North America, Europe, Asia, the east coast of South America and Africa, and southern Australia/Oceania (Fig. S7d). However, FBC and BBC in the top 1 m widely increased due to the positive effects of temperature and substrate availability (Figs. 6 and 7c–f). Therefore, negative correlations of FBC and BBC with SM in the top 1 m at middle and low latitudes in North America, Europe, Asia, the east coast of South America and Africa, and southern Australia/Oceania resulted from the facilitating effects of temperature and substrate availability on FBC and BBC in the top 1 m despite decreasing water availability.

4.4 Future improvements

Although the CLM-Microbe model can well reproduce the global distribution of C in vegetation, soil, and microbes, four key improvements are identified for future work. First, soil and microbial processes along soil profiles need to be better defined. Soil and microbial variables such as DOC, SOC, FBC, and BBC in 0–30 cm were better simulated than those in 0–1 m (Table 1; Figs. 1–2), indicating that soil and microbial processes in the deeper soil profile are not adequately modeled. Therefore, better defining soil and microbial processes with depth can help improve the model efficacy in capturing soil and microbial processes and further reduce uncertainties in future projections of the C cycle. Second, land-use change needs to be considered in future work. In addition to changes in environmental factors, land-use change also has profound influences on the plant, soil, and microbial processes. Drastic changes in vegetation, soil, and microbial processes due to land-use change can occur at small scales, and the spatial pattern of those processes can also be changed (Pascual et al., 1997; Sampaio et al., 2007; Stevenson et al., 2016). Therefore, considering the impacts of land-use change in the CLM-Microbe model can help improve the model efficiency in capturing spatial patterns of C density and stocks in terrestrial ecosystems. The global biogeographic patterns of soil microbes and their functions have been recognized (Xu et al., 2020). This modeling study has made progress toward a full investigation of microbial patterns and mechanisms, and a community-wide microbial data system is needed to facilitate more data–model integration to improve microbial models. Lastly, factorial analysis to attribute the variations in terrestrial C fluxes will be addressed in our future work. Variations in terrestrial C fluxes and pools are driven by multiple environmental change factors that contribute individually or in combination. Attributing the variations in terrestrial C fluxes and pools to environmental change factors is important for the understanding of terrestrial C flux and pool dynamics (Xu et al., 2010).

5 Conclusions

The ESMs incorporating microbial processes are expected to represent uncertainties in the terrestrial C cycle more comprehensively. The CLM-Microbe model can reproduce the distribution of vegetation (GPP, NPP, and VegC), soil (HR, SR, DOC, and SOC), and microbial (FBC, BBC, and MBC) variables. In addition, microbial fluxes (HR and SR) and pools (FBC and BBC in top 30 cm and 1 m) increased from 1901 to 2016. We observed increases of such variables in most of the land but slight decreases of FBC and BBC in the top 1 m in Australia/Oceania. The increases in HR, SR, and fungal and bacterial biomass C were closely associated with increasing vegetation C input and SM and ST in the top 1 m.

This study represents one of the first attempts to simulate the spatial and temporal variations in C fluxes and pool sizes of soil microbes during the last century using a microbial-explicit model – the CLM-Microbe model. As the community is moving towards a microbial-explicit Earth system model, this study provides robust support for microbial model development and application for predicting microbial roles in the C–climate feedback. The variations in soil microbial community over historical periods and across space simulated by the CLM-Microbe model provide a crucial foundation to study the impacts of soil microbes on terrestrial biogeochemical processes.