Global soil consumption of atmospheric carbon monoxide: an analysis using a process-based biogeochemistry model

Carbon monoxide (CO) plays an important role in controlling the oxidizing capacity of the atmosphere by reacting with OH radicals that affect atmospheric methane (CH₄) dynamics. We develop a process-based biogeochemistry model to quantify the CO exchange between soils and the atmosphere with a 5 min internal time step at the global scale. The model is parameterized using the CO flux data from the field and laboratory experiments for 11 representative ecosystem types. The model is then extrapolated to global terrestrial ecosystems using monthly climate forcing data. Global soil gross consumption, gross production, and net flux of the atmospheric CO are estimated to be from −197 to −180, 34 to 36, and −163 to −145 Tg CO yr⁻¹ (1 Tg = 10¹² g), respectively, when the model is driven with satellite-based atmospheric CO concentration data during 2000–2013. Tropical evergreen forest, savanna and deciduous forest areas are the largest sinks at 123 Tg CO yr⁻¹. The soil CO gross consumption is sensitive to air temperature and atmospheric CO concentration, while the gross production is sensitive to soil organic carbon (SOC) stock and air temperature. By assuming that the spatially distributed atmospheric CO concentrations (∼ 128 ppbv) are not changing over time, the global mean CO net deposition velocity is estimated to be 0.16–0.19 mm s⁻¹ during the 20th century. Under the future climate scenarios, the CO deposition velocity will increase at a rate of 0.0002–0.0013 mm s⁻¹ yr⁻¹ during 2014–2100, reaching 0.20–0.30 mm s⁻¹ by the end of the 21st century, primarily due to the increasing temperature. Areas near the Equator, the eastern US, Europe and eastern Asia will be the largest sinks due to optimum soil moisture and high temperature. The annual global soil net flux of atmospheric CO is primarily controlled by air temperature, soil temperature, SOC and atmospheric CO concentrations, while its monthly variation is mainly determined by air temperature, precipitation, soil temperature and soil moisture.

1 Introduction

Carbon monoxide (CO) plays an important role in controlling the oxidizing capacity of the atmosphere by reacting with OH radicals (Logan et al., 1981; Crutzen, 1987; Khalil and Rasmussen, 1990; Prather et al., 1995; Prather and Ehhalt, 2001). CO in the atmosphere can directly and indirectly influence the fate of critical greenhouse gases such as methane (CH₄) and ozone (O₃) (Tan and Zhuang, 2012). Although CO itself absorbs only a limited amount of infrared radiation from the Earth, the cumulative indirect radiative forcing of CO may be even larger than that of the third powerful greenhouse gas, nitrous oxide (N₂O, Myhre et al., 2013). Current estimates of global CO emissions from both anthropogenic and natural sources range from 1550 to 2900 Tg CO yr⁻¹, which are mainly from anthropogenic and natural direct emissions and from the oxidation of methane and other volatile organic compounds (VOCs) (Prather et al., 1995; Khalil et al., 1999; Bergamaschi et al., 2000; Prather and Ehhalt, 2001; Stein et al., 2014). Chemical consumption of CO by atmospheric OH and the biological consumption of CO by soil microbes are two major sinks of the atmospheric CO (Conrad, 1988; Lu and Khalil, 1993; Yonemura et al., 2000; Whalen and Reeburgh, 2001).

Soils are globally considered as a major sink for CO due to microbial activities (Whalen and Reeburgh, 2001; King and Weber, 2007). A diverse group of soil microbes including carboxydotrophs, methanotrophs and nitrifiers are capable of oxidizing CO (King and Weber, 2007). Annually, 10–25 % of total earth surface CO emissions were consumed by soils (Sanhueza et al., 1998; King, 1999a; Chan and Steudler, 2006). Potter et al. (1996) reported the global soil consumption to be from −50 to −16 Tg CO yr⁻¹ (negative values represent the uptake from the atmosphere to soil), by using a single-box model over the upper 5 cm of soils. All existing estimates have large uncertainties and range from −640 to −16 Tg CO yr⁻¹ (Sanhueza et al., 1998; King, 1999b; Bergamaschi et al., 2000). Similarly, the estimates of CO dry deposition velocities also have large uncertainties and range from 0 to 4.0 mm s⁻¹; here, positive values represent deposition to soils (King, 1999a; Castellanos et al., 2011). Soils also produce CO mainly via abiotic processes such as thermal degradation and photo-degradation of organic matter or plant materials (Conrad and Seiler, 1985; Tarr et al., 1995; Schade and Crutzen, 1999; Derendorp et al., 2011; Lee et al., 2012; van Asperen et al., 2015; Fraser et al., 2015; Pihlatie et al., 2016), except for a few cases of anaerobic formation. Photo-degradation is identified as radiation-dependent degradation due to absorbing radiation (King et al., 2012). Thermal degradation is identified as the temperature-dependent degradation of carbon in the absence of radiation and possibly oxygen (Derendorp et al., 2011; Lee et al., 2012; van Asperen et al., 2015; Pihlatie et al., 2016). These major soil CO production processes, together with soil CO consumption processes, have not been adequately considered in global soil CO budget estimates.

To date, most top–down atmospheric models have applied a dry deposition scheme based on the resistance model of Wesely (1989). Such schemes provided a wide range of dry deposition velocities (Stevenson et al., 2006). Only a few models (MOZART-4, Emmons et al., 2010; CAM-chem, Lamarque et al., 2012) have extended their dry deposition schemes with a parameterization for CO and H₂ uptake through oxidation by soil microbes, following the work of Sanderson et al. (2003), which was based on extensive measurements from Yonemura et al. (2000). Potter et al. (1996) developed a bottom–up model to simulate CO consumption and production at the global scale. Their model is a single-box model that only considers the top 5 cm depth of soil and does not have explicit microbial factors, and therefore might have underestimated CO consumption (Potter et al., 1996; King, 1999a). Current bottom–up CO modeling approaches are mostly based on a limited number of CO in situ observations or laboratory studies to quantify regional and global soil consumption (Potter et al., 1996; Sanhueza et al., 1998; Khalil et al., 1999; King, 1999a; Bergamaschi et al., 2000; Prather and Ehhalt, 2001). To our knowledge, no detailed process-based model of soil–atmospheric exchange of CO has been published in the last 15 years. One reason is that there is an incomplete understanding of biological processes of uptake (King and Weber, 2007; Vreman et al., 2011; He and He, 2014; Pihlatie et al., 2016). Another reason is that there is a lack of long-term CO flux measurements for different ecosystem types to calibrate and evaluate the models. CO flux measurements are mostly from short-term field observations or laboratory experiments (e.g., Conrad and Seiler, 1985; Funk et al., 1994; Tarr et al., 1995; Zepp et al., 1997; Kuhlbusch et al., 1998; Moxley and Smith, 1998; Schade et al., 1999; King and Crosby, 2002; Varella et al., 2004; Lee et al., 2012; Bruhn et al., 2013; van Asperen et al., 2015). The first study to report long-term and continuous field measurements of CO flux over grasslands using a micrometeorological eddy covariance (EC) method is Pihlatie et al. (2016).

To improve the quantification of the global soil CO budget for the period 2000–2013 and the CO deposition velocity for the 20th and 21st centuries, this study developed a CO dynamics module (CODM) embedded in a process-based biogeochemistry model, the Terrestrial Ecosystem Model (TEM) (Zhuang et al., 2003, 2004, 2007). The CODM was then calibrated and evaluated using laboratory experiments and field measurements for different ecosystem types. The atmospheric CO concentration data from MOPITT (Gille, 2013) were used to drive model simulations from 2000 to 2013. A set of century-long simulations of 1901–2100 were also conducted using the atmospheric CO concentrations estimated with an empirical function (Badr and Probert, 1994; Potter et al., 1996). Finally, the effects of multiple forcings on the global CO consumption and production, including the changes in climate and atmospheric CO concentrations at the global scale, were evaluated with the model.

2 Method

2.1 Overview

We first developed a soil CO dynamics module (CODM) on a daily time step that considers (1) the soil–atmosphere CO exchange and diffusion process between soil layers, (2) the consumption by soil microbial oxidation, (3) the production by soil chemical oxidation, and (4) the effects of temperature, soil moisture, soil CO substrate and surface atmospheric CO concentration on these processes. Second, we used the observed soil temperature and moisture to evaluate the TEM hydrology module and the soil thermal module in order to estimate soil physical variables. Then we used the data from laboratory experiments and CO flux measurements to parameterize the model using the Shuffled Complex Evolution (SCE-UA) method (Duan et al., 1993). Finally, the model was extrapolated to the globe at a 0.5^∘ by 0.5^∘ resolution. We conducted three sets of model experiments to investigate the impact of climate and atmospheric CO concentrations on soil CO dynamics: (1) simulations for 2000–2013 with MOPITT satellite atmospheric CO concentration data; (2) simulations for 1901–2100 with constant atmospheric CO concentrations estimated from an empirical function and the historical climate data (1901–2013) and three future climate scenarios (2014–2100); and (3) eight sensitivity simulations by increasing and decreasing (a) constant CO surface concentrations by 30 %, (b) SOC by 5 %, (c) precipitation by 20 % and (d) air temperature by 3^∘C for each pixel, respectively, while holding other forcing data as they were, during 1999–2000.

Figure 1 The model framework includes a carbon and nitrogen dynamics module (CNDM), a soil thermal module (STM) from Terrestrial Ecosystem Model (TEM) 5.0 (Zhuang et al., 2001, 2003), a hydrological module (HM) based on a land surface module (Bonan, 1996; Zhuang et al., 2004), and a carbon monoxide dynamics module (CODM). The detailed structure of CODM includes land surface CO concentration as the top boundary and 30 1 cm thick layers (in total, 30 cm) where consumption and production take place.

2.2 Carbon monoxide dynamics module (CODM)

Embedded in the TEM (Fig. 1), the CODM is mainly driven by (1) soil organic carbon availability based on a carbon and nitrogen dynamics module (CNDM) (Zhuang et al., 2003); (2) a soil temperature profile from a soil thermal module (STM) (Zhuang et al., 2001, 2003); and (3) a soil moisture profile from a hydrological module (HM) (Bonan, 1996; Zhuang et al., 2004). The net exchange of CO between the atmosphere and soil is determined by the mass balance approach (net flux = total production – total oxidation – total soil CO concentration change). According to previous studies, we separated active soils (top 30 cm) for CO consumption and production into 1 cm thickness layers (King, 1999a, b; Whalen and Reeburgh, 2001; Chan and Steudler, 2006). Between the soil layers, the changes in CO concentrations were calculated as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \left(C\right(t,i\left)\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D\left(t,i\right){\displaystyle \frac{\partial \left(C\left(t,i\right)\right)}{\partial z}}\right)+P\left(t,i\right)-O\left(t,i\right),\end{array}$$

where C(ti) is the CO concentration (mg m⁻³) in layer i and at time t. z is the depth of the soil (m). D(t,i) is the diffusion coefficient (m² s⁻¹) for layer i. P(t,i) is the CO production rate (mg m⁻³ s⁻¹) and O(t,i) is the CO consumption rate (mg m⁻³ s⁻¹). D(t,i) is calculated using the method from Potter et al. (1996), which is a function of soil temperature, soil texture and soil moisture. The upper boundary condition is the atmospheric CO concentration, which is estimated with an empirical function of latitude (Potter et al., 1996) or directly measured by the MOPITT satellite during 2000–2013. The lower boundary condition is assumed to have no diffusion exchange with the layer underneath. This partial differential equation (PDE) is solved using the Crank–Nicolson method for a less time-step-sensitive solution.

The CO consumption was modeled in unsaturated soil pores as

$$\begin{array}{}\text{(2)}& {\displaystyle }O\left(t,i\right)=V_{\text{max}}\cdot f_{\mathrm{1}}\left(C\left(t,i\right)\right)\cdot f_{\mathrm{2}}\left(T\left(t,i\right)\right)\cdot f_{\mathrm{3}}\left(M\left(t,i\right)\right),\end{array}$$

where V_max is the ecosystem-specific maximum oxidation rate and was estimated previously ranging from 0.3 to 11.1 µg CO g⁻¹ h⁻¹ for different ecosystems (Whalen and Reeburgh, 2001). f_i represents the effects of soil CO concentration C(t,i), temperature T(t,i) and moisture M(t,i) on the CO soil consumption. Considering the CO consumption as the result of microbial activities, we calculated f₁(C(t,i)), f₂(T(t,i)) and f₃(M(t,i)) in a similar way to Zhuang et al. (2004):

$$\begin{array}{ll}\text{(2a)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{1}}\left(C\left(t,i\right)\right)={\displaystyle \frac{C\left(t,i\right)}{C\left(t,i\right)+k_{\mathrm{CO}}}},\text{(2b)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{2}}\left(T\left(t,i\right)\right)=Q_{\mathrm{10}}^{\frac{T\left(t,i\right)-T_{\text{ref}}}{\mathrm{10}}},{\displaystyle }& {\displaystyle }{f}_{\mathrm{3}}\left(M\left(t,i\right)\right)=\\ \text{(2c)}& {\displaystyle }& {\displaystyle \frac{(M\left(t,i\right)-M_{\text{min}})(M\left(t,i\right)-M_{\text{max}})}{\left(M\left(t,i\right)-M_{\text{min}}\right)\left(M\left(t,i\right)-M_{\text{max}}\right)-(M\left(t,i\right)-M_{\text{opt}}{)}^{\mathrm{2}}}},\end{array}$$

where f₁(C(t,i)) is a multiplier that enhances the oxidation rate with increasing soil CO concentrations using a Michaelis–Menten function with a half-saturation constant k_CO, and their values were previously estimated ranging from 5 to 51 µl CO l⁻¹ for different ecosystems (Whalen and Reeburgh, 2001); f₂(T(t,i)) is a multiplier that enhances the CO oxidation rate with increasing soil temperature using a Q₁₀ function with Q₁₀ coefficients (Whalen and Reeburgh, 2001). T_ref is the reference temperature; units are ^∘C (Zhuang et al., 2004, 2013). f₃(M(t,i)) is a multiplier to estimate the biological limiting effect that diminishes the CO oxidation rates if the soil moisture is not at an optimum level (M_opt). M_min, M_max and M_opt are the minimum, maximum and optimum volumetric soil moisture of oxidation reaction, respectively. Equation (2b) will overestimate the CO consumption at high temperature because in reality the CO consumption will decrease when temperature is higher than optimum temperature, while f₂ will keep increasing with rising temperature. However, the CO consumption is constrained by the CO production, and Eq. (1) is used to represent this constraint.

We modeled the CO production rate (P(t,i)) as a process of chemical oxidation constrained by the soil organic carbon (SOC) decay (Conrad and Seiler, 1985; Potter et al., 1996; Jobbagy and Jackson, 2000; van Asperen et al., 2015):

$$\begin{array}{}\text{(3)}& {\displaystyle }P\left(t,i\right)=P_{\text{r}}\left(t,i\right)\cdot E_{\text{SOC}}\cdot C_{\text{SOC}}\left(t\right)\cdot F_{\text{SOC}},\end{array}$$

where P_r(t,i) is a reference soil CO production rate which has been normalized to the rate at reference temperature (the production rate at temperature (t,i) divided by the production rate at the reference temperature), which is affected by soil moisture and soil temperature (Conrad and Seiler, 1985; van Asperen et al., 2015). E_SOC is an estimated nominal CO production factor of 3.5 ± 0.9 × 10⁻⁹ mg CO m⁻² s⁻¹ per g SOC m⁻² (to 30 cm soil depth) (Potter et al., 1996). C_SOC(t) is a SOC content (mg m⁻²), which is provided by a CNDM in the TEM. F_SOC is a constant fraction of top 30 cm SOC compared to the total amount of SOC, which is 0.33 for shrubland areas, 0.42 for grassland areas and 0.50 for forest areas, respectively (Jobbagy and Jackson, 2000). P_r(t,i) was calculated as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{P}_{\text{r}}\left(t,i\right)=\exp \left(f_{\mathrm{4}}\left(M\left(t,i\right)\right)\cdot {\text{Ea}}_{\text{ref}}/R\cdot \right.\\ \text{(3a)}& {\displaystyle }& {\displaystyle }\left.\left({\displaystyle \frac{\mathrm{1}}{\mathrm{273.15}+{\text{PT}}_{\text{ref}}}}-{\displaystyle \frac{\mathrm{1}}{T\left(t,i\right)+\mathrm{273.15}}}\right)\right),\text{(3b)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{4}}\left(M\left(t,i\right)\right)={\displaystyle \frac{{\text{PM}}_{\text{ref}}}{M\left(t,i\right)+{\text{PM}}_{\text{ref}}}},\end{array}$$

where Eq. (3a) is derived from the Arrhenius equation for chemical reactions and normalized using the reference temperature PT_ref. Ea_ref∕R is the reference activation energy divided by gas constant R; units are K. f₄(M(t,i)) is the multiplier that reduces activation energy using a regression approach based on the laboratory experiment of moisture influences on CO production (Conrad and Seiler, 1985). PM_ref is the reference volumetric soil moisture, ranging from 0.01 to 0.5 volume volume⁻¹ (v v⁻¹). We assumed the thermal degradation to be the main CO producing process due to a lack of photo-degradation data and it being hard to distinguish photo-degradation from observations. In order to reduce the bias from the thermal degradation to the total abiotic degradation, Eq. (3a) is parameterized by comparing with the total production rate. For instance, P_r(t,i) calculation can perfectly fit the experiment results in Van Asperen et al. (2015) with proper PT_ref (18 ^∘C), Ea_ref∕R (14 000 K), and PM_ref (0.5 v v⁻¹).

Table 1 Model parameterization sites for the thermal and hydrology modules (site nos. 1–4) and for the CODM (site nos. 5–13).

The CO deposition velocity was modeled in the same way as Eq. (19.1) in Seinfeld and Pandis (1998):

$$\begin{array}{}\text{(4)}& v_{\text{d}}=-F_{\text{net}}/C_{\text{CO,air}},\end{array}$$

where v_d is the CO deposition velocity (mm s⁻¹). F_net is the model estimated CO net flux rate (mg CO m⁻² day⁻¹). C_CO,air is the CO surface concentration (ppbv). C_CO,air can be MOPITT CO surface concentration data or the derived CO surface concentrations using the same method as Potter et al. (1996). Positive values of v_d represent soil uptake (deposition from air to soils) and negative values represent soil emissions.

Figure 2 Evaluation of thermal and hydrology module at four sites: (a) Boreal evergreen needleleaf forests; (b) temperate deciduous broadleaf forests. (1) shows the soil temperature comparison between the model simulations (gray line) and observations (black line) and (2) shows the soil moisture comparison between the model simulations (gray line) and observations (black line). Specifically, the volumetric soil moisture is converted from the water content reflectometry (WCR) probe output period using an empirical calibration function of Bourgeau-Chavez et al. (2012) for the 5–30 cm layer. Some of them resulted in calculations of values greater than 100 % volumetric soil moisture (VSM) in the Nakai et al. (2013) study. Our model estimated high VSM (close to 80 %) is due to top 10 cm moss in the model, which has a saturation VSM of 0.8. Evaluation of the thermal and hydrology module at four sites: (c) tropical moist forest, (d) tropical forest–savanna. (1) shows the soil temperature comparison between the model simulations (gray line) and observations (black line) and (2) shows the soil moisture comparison between the model simulations (gray line) and observations (black line).

2.3 Model parameterization and extrapolation

The model parameterization was conducted in two steps: (1) thermal and hydrology modules embedded in the TEM were revised, calibrated and evaluated by running a model driven by corresponding local meteorological or climatic data at four representative sites, including boreal forest, temperate forest, tropical forest and savanna (Table 1, site nos. 1 to 4, Fig. 2) to minimize model–data mismatch in terms of soil temperature and moisture. (2) The CODM module was parameterized by running the TEM for observational periods driven with the corresponding local meteorological or climatic data at each reference site (Table 1, Fig. 3), and using the Shuffled Complex Evolution Approach in the R language (SCE-UA-R) (Duan et al., 1993) to minimize the difference between the simulated and observed net CO flux. Eleven parameters including k_CO, V_max, T_ref, Q₁₀, M_min, M_max, M_opt, E_SOC, Ea_ref∕R, PM_ref and PT_ref were optimized (Table 2). F_SOC was not involved in the calibration process. Parameter priors were decided based on previous studies (Conrad and Seiler, 1985; King, 1999b; Whalen and Reeburgh, 2001; Zhuang et al., 2004). The SCE-UA-R was used for site nos. 6, 8, 10, 11 (Table 1). In parameter ensemble simulations, we have run 50 times SCE-UA-R with 10 000 maximum loops for each site, and all of them have reached a stable state before the end of the loops. For wetlands, the only available data for calibration are from site no. 12. We used the trial-and-error method to make the simulated results in the range of observed flux rates, with a 10 % tolerance. For tropical sites, since the tropical savanna vegetation type is treated as a combination type of tropical forest and grassland in our simulations, we first used site no. 13 to set priors to fit the experiment results with a 10 % tolerance and then evaluated it by running our model and comparing it with the site no. 7 results. Site nos. 9 and 5 were used to evaluate our model results for temperate forest and grassland. Besides the observed climatic and soil property data, we used the ERA-Interim reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011), AmeriFlux observed meteorology data (US-Dk3, Novick et al., 2016; CA-Man, Amiro, 2016) and reanalysis climatic data from the Climatic Research Unit (CRU, Harris et al., 2013) to fill the missing environmental data. To sum up, parameters for various ecosystem types in Table 2 were the final results of our parameterization. Model parameterization was conducted for ecosystem types including boreal forest, temperate coniferous forest, temperate deciduous forest, and grassland using SCE-UA-R. In contrast, for tropical forest and wet tundra, we used a trial-and-error method to adjust parameters to allow the model simulation to best fit the observed data. Due to limited data availability, we assumed temperate evergreen broadleaf forests have the same parameters as temperate deciduous forest ecosystems.

Table 2 Ecosystem-specific parameters in the CODM^a.

^a k_CO is the half-saturation constant for soil CO concentration; V_max is the specific maximum CO oxidation rate; T_ref is the reference temperature to account for the soil temperature effects on the CO consumption; Q10 is the ecosystem-specific Q10 coefficient to account for soil temperature effects on the CO consumption; M_min, M_max, and M_opt are the minimum, optimum, and maximum volumetric soil moisture of oxidation reaction to account for soil moisture effects on the CO consumption; E_SOC is an estimated nominal CO production factor, similar to Potter et al. (1996) (10⁻⁴ mg CO m⁻² d⁻¹ per g SOC m⁻²); F_SOC is a constant fraction of top 20 cm SOC compared to the total amount of SOC to account for SOC effects on the CO production; Ea_ref∕R is the ecosystem-specific activation energy divided by a gas constant to account for the reaction rate of production; PM_ref is the reference moisture to account for soil temperature effects on the CO production; PT_ref is the reference temperature to account for soil temperature effects on the CO production.

Figure 3 Parameter ensemble experiment results: each parameter has 50 calibrated values generated from running SCE-UA-R 50 times independently. Parameters are normalized to their largest potential values described in Table 2. (a1) and (a2) are temperate coniferous forest normalized parameter distribution boxplots and CO flux comparisons between the model simulations (solid line, using the mean value of parameters) and observations (green diamond; red lines represent the error bar, site no. 8), respectively. For each box, line top, box top, horizontal line inside box, box bottom and line bottom represent maximum, third quartile, median, first quartile and minimum of 50 parameter values. Red dot represents the mean value of 50 parameter values. (b1) and (b2) are plots for temperate deciduous forest (site no. 11). Parameter ensemble experiment results: each parameter has 50 calibrated values generated from running SCE-UA-R 50 times independently. Parameters are normalized to their largest potential values described in Table 2. (c1) and (c2) are boreal forest normalized parameter distribution boxplots and CO flux comparisons between the model simulations (solid line, using the mean value of parameters) and observations (green diamond; red lines represent the error bar, site no. 12), respectively. For each box, line top, box top, horizontal line inside box, box bottom and line bottom represent maximum, third quartile, median, first quartile and minimum of 50 parameter values. Red dot represents the mean value of 50 parameter values. (d1) and (d2) are for grassland (site no. 6). Grassland observation data are the sum of hourly observations, so the error bar represents the standard deviation.

2.4 Data organization

To get the spatially and temporally explicit estimates of the CO consumption, production and net flux at the global scale, we used the data of land cover, soils, climate and leaf area index (LAI) from various sources at a spatial resolution of 0.5^∘ latitude × 0.5^∘ longitude to drive the TEM. The land cover data include the potential vegetation distribution (Melillo et al., 1993) and soil texture (Zhuang et al., 2003), which were used to assign vegetation- and texture-specific parameters to each grid cell.

Figure 4 Historical global land surface (excluding the Antarctic area and ocean area) mean climate, and simulated global mean soil moisture, soil temperature and SOC for the period 1901–2013.

For the simulation of the period 1901–2013, the monthly air temperature, precipitation, cloud fraction and vapor pressure data sets from CRU were used to estimate the soil temperature, soil moisture and SOC with the TEM (Fig. 4). The monthly LAI data from the TEM were required to simulate soil moisture (Zhuang et al., 2004). During this time period, we used an empirical function of latitude which was derived from the observed latitudinal distribution of tropospheric carbon monoxide (Badr and Probert, 1994) to calculate the static CO surface concentration distribution (Eq. 7, Potter et al., 1996):

$$\begin{array}{ll}{\displaystyle }{C}_{\text{CO,air}}& {\displaystyle }=\mathrm{82.267856}+\mathrm{0.8441503}L+\mathrm{1.55934}\times {\mathrm{10}}^{-\mathrm{2}}{L}^{\mathrm{2}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }+\mathrm{2.37}\times {\mathrm{10}}^{-\mathrm{5}}{L}^{\mathrm{3}}-\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{6}}{L}^{\mathrm{4}},\end{array}$$

Figure 5 CO surface concentration data from the MOPITT satellite (ppbv): (a) the global mean CO surface concentrations from MOPITT during 2000–2013; (b) the CO annual surface concentrations from both MOPITT and empirical functions (Potter et al., 1996).

where C_CO,air is the derived surface CO concentration (ppbv), and L represents latitudes with negative degrees for the Southern Hemisphere and positive degrees for the Northern Hemisphere. We also used the atmospheric CO data from the MOPITT satellite during 2000–2013 (Fig. 5). We averaged day-time and night-time monthly mean values of CO surface level 3 retrieval data (variables mapped on 0.5^∘ latitude × 0.5^∘ longitude grid scales with a monthly time step, Gille, 2013) to represent the CO surface concentration level in each month. The pixels with missing values were filled with the average values of those pixels that were inside 1.5 times of the distance between the missing-value pixel and the nearest pixel with values. These global mean values shown in Fig. 5 do not include ocean surfaces; thus, there are differences between our surface CO concentration results and Yoon and Pozzer's (2014) report in 2014, which is as low as 99.8 ppb. From 2014 to 2100, we used the Intergovernmental Panel on Climate Change (IPCC) future climate scenarios from Representative Concentration Pathways (RCPs) climate forcing data sets RCP2.6, RCP4.5 and RCP8.5 (Fig. 6). RCP2.6, 4.5 and 8.5 data sets are future climate projections with anthropogenic greenhouse gas emission radiative forcings of 2.6, 4.5 and 8.5 W m⁻², respectively, by 2100. Since RCPs did not have water vapor pressure data, we used the specific humidity and sea level air pressure from the RCPs and elevation of surface to estimate the monthly surface vapor pressure (Seinfeld and Pandis, 1998).

Figure 6 Global land surface (excluding the Antarctic area and ocean area) mean climate from the RCP2.6, RCP4.5 and RCP8.5 data sets and simulated mean soil temperature, moisture and SOC: (a–g) are land surface air temperature (^∘C), soil temperature (^∘C), precipitation (mm), soil moisture (%), surface water vapor pressure (hpa), cloud fraction (%), and SOC (mg m⁻²), respectively.

2.5 Model experiment design

We conducted two sets of core simulations and eight sensitivity test simulations for a historical period. The two core sets of simulations were driven with the MOPITT CO surface concentration data for the period 2000–2013 (experiment E1) and with the spatially distributed CO surface concentrations assuming them as constant over time estimated from an empirical function of latitude for the period 1901–2100 (experiment E2), respectively. Specifically, in experiment E2 we used the CRU climate forcing for the historical period 1901–2013 and the climate data of RCP2.6, RCP4.5 and RCP8.5 for different future scenarios to examine the responses of CO flux to changing climates. Eight sensitivity simulations were driven with varying different forcing variables while keeping others as they were: (1) with constant CO surface concentrations ± 30 %, (2) SOC ± 5 %, (3) precipitation ± 20 % and (4) air temperature ± 3 ^∘C for each pixel, respectively, during 1999–2000 (E3).

3 Results

3.1 Site evaluation

Both the magnitude and variation of the simulated soil temperature and moisture from cold areas to warm areas compare well to the observations (Fig. 2). The magnitude of the simulated CO flux is comparable and correlated with the observations (R is about 0.5, p-value < 0.001, Fig. 3a2, b2, c2, d2). Estimated CO fluxes for different ecosystem types range from −28.4 to 1.7 mg CO m⁻² day⁻¹, and the root mean square error (RMSE) between simulation and observation at all sites is below 1.5 mg CO m⁻² day⁻¹. RMSE for site no. 7 is bigger than 2.0 mg CO m⁻² day⁻¹ when compared with transparent chamber observations. For the boreal forest site, we only had eight acceptable points in 1994 and 1996 (Fig. 3c2).

Figure 7 Global mean soil CO consumption, production and net flux: (a) the annual time series during 2000–2013 and (b) the latitudinal distribution during 2000–2013.

Table 3 Regional soil CO consumption, net flux and production (Tg CO yr⁻¹) during 2000–2013.

3.2 Global soil CO dynamics during 2000–2013

Using the MOPITT CO surface concentration data during 2000–2013 (E1), the estimated mean soil CO consumption, production and net flux (positive values indicate CO emissions from soils to the atmosphere) are from −197 to −180, 34 to 36 and −163 to −145 Tg CO yr⁻¹, respectively (Fig. 7a). The consumption is about 4 times larger than the production. The annual consumption and net flux trends follow the atmospheric CO concentration trends (Figs. 5b, 7a), with a small interannual variability (< 10 %). The latitudinal distributions of the consumption, production and net fluxes share the same spatial pattern. Around 20^∘ S–20^∘ N and 20–60^∘ N are the largest and second largest areas for production and consumption, while the 45^∘ S–45^∘ N area accounts for nearly 90 % of the total consumption and production (Fig. 7b, Table 3). The Southern Hemisphere and Northern Hemisphere have 41 and 59 % of the total consumption, and 47 and 53 % of the total production, respectively (Table 3). The highest rates of the consumption and production are located in areas close to the Equator, and the consumption from areas such as the eastern US, Europe and eastern Asia is also high (< −1000 mg m⁻² yr⁻¹) (Fig. 8a, b). Global soils serve as an atmospheric CO sink (Fig. 8c). Some areas, such as the western US and southern Australia, are CO sources, all of which are grasslands or experiencing dry climate. Tropical evergreen forests are the largest sinks, consuming 86 Tg CO yr⁻¹, and tropical savanna and deciduous forest are the second and third largest sinks, consuming a total of 37 Tg CO yr⁻¹ (Table 4). These three ecosystems account for 66 % of the total consumption. Tropical evergreen forests are also the largest source of soil CO production, producing 16 Tg CO yr⁻¹, while tropical savanna has a considerable production of 6 Tg CO yr⁻¹ (Table 4). Moreover, tropical areas, including forested wetlands, forested floodplain and evergreen forests, are most efficient for the CO consumption, ranging from −18 to −13 mg CO m⁻² day⁻¹. They are also most efficient for the CO production, at over 2 mg CO m⁻² day⁻¹ (Table 4, calculated by fluxes divided by area).

Table 4 Annual total soil CO consumption, net flux and production in different ecosystems during 2000–2013 (E1) and mean CO deposition velocity in different ecosystems during 1901–2013 (E2).

Figure 8 Global annual mean soil CO fluxes (mg CO m⁻² yr⁻¹) during 2000–2013 using the MOPITT CO atmospheric surface concentration data.

Figure 9 Global mean annual time series of CO deposition velocity (mm s⁻¹) using constant in time and spatially distributed CO concentration data during 1901–2013 (a) and under the future climate scenarios RCP2.6, RCP4.5 and RCP8.5 during 2014–2100 (b).

3.3 Global soil CO dynamics during 1901–2100

Using the constant CO surface concentration, the estimated global mean CO deposition velocities are 0.16–0.19 mm s⁻¹ for the period 1901–2013. For the period 2014–2100, the deposition velocities are 0.18–0.21, 0.18–0.24 and 0.17–0.31 for the RCP2.6, 4.5 and 8.5 scenarios, respectively (Fig. 9). During 2014–2100, there are significant trends of increasing deposition velocities for nearly all scenarios (Fig. 9). The rates of increase are 0.0002, 0.0005 and 0.0013 mm s⁻¹ yr⁻¹, and will reach 0.20, 0.23 and 0.30 mm s⁻¹ by the end of the 21st century for the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively (Fig. 9). These increasing trends are similar to the air temperature increasing trends (Fig. 6a). The global distribution patterns of the CO deposition velocity are similar to the net flux distribution for the period 2000–2013, but there are significant differences between the 1901–2013, RCP2.6, RCP4.5 and RCP8.5 scenarios (Fig. 10). The deposition velocities increase from the RCP2.6 to RCP8.5 and larger than that in the historical periods in areas near the Equator (Fig. 10). Areas near the Equator and eastern Asia become big sinks of the atmospheric CO, while the northeastern US becomes a small source in the 21st century (Fig. 10). Different vegetation types have a large range of the deposition velocity, from 0.008 to 1.154 mm s⁻¹ (Table 4). The tropical forested wetland, tropical forested floodplain and tropical evergreen forest have the top three largest deposition velocities of 1.154, 1.117 and 0.879 mm s⁻¹, respectively, while desert, short grasslands, and wet tundra have the smallest deposition velocities of 0.008, 0.010 and 0.015 mm s⁻¹, respectively.

Figure 10 Global annual mean CO deposition velocity using constant in time and spatially distributed CO concentration data (mm s⁻¹) (a) during 1901–2013 and (b), (c), (d) under the future climate scenarios RCP2.6, RCP4.5 and RCP8.5 during 2014–2100, respectively.

3.4 Sensitivity test

The soil CO consumption is most sensitive (changing 29 %) to air temperature, while the production is most sensitive to both air temperature (changing up to 36 %) and SOC (5 %). The net CO fluxes have similar sensitivities to the consumption. The annual CO consumption, production and net flux follow the change in air temperature (Table 5). In addition, a 30 % change in precipitation will not lead to large changes in the CO flux (< 3 %).

Table 5 Sensitivity of the global CO consumption, net flux and production (Tg CO yr⁻¹) to the changes in atmospheric CO, soil organic carbon (SOC), precipitation (Prec) and air temperature (AT).

4 Discussion

4.1 Comparison with other studies

Previous studies estimated a wide range of the global CO consumption from −16 to −640 Tg CO yr⁻¹. Our estimates are from −197 to −180 Tg CO yr⁻¹ for 2000–2013 using the MOPITT satellite CO surface concentration data. Previous studies also provided a large range for the CO production from 0 to 7.6 mg m⁻² day⁻¹ (reviewed in Pihlatie et al., 2016). Our results showed the averaged CO production ranging from 0.01 to 2.29 mg m⁻² day⁻¹. The existing estimates of the CO deposition velocities for different vegetation types ranged from 0.0 to 4.0 mm s⁻¹, while our simulations showed an averaged CO deposition velocity ranging from 0.006 to 1.154 mm s⁻¹ for different vegetation types. The large uncertainty of these estimates is mainly due to different considerations of the microbial activities, the depth of the soil, and the parameters in the model. In contrast to the estimates of −57 to −16 Tg CO yr⁻¹ which were based on the top 5 cm of soils (Potter et al., 1996), our estimates considered 30 cm soils, which were used in Whalen and Reeburgh (2001). In addition, we used a thinner layer division (1 cm for each layer) for the diffusion process, and used the Crank–Nicolson method to solve partial differential equations to avoid time step influences. We also included the microbial CO oxidation process to remove the CO from the soil and considered the effects of soil moisture, soil temperature, vegetation type and soil CO substrate on microbial activities. Our soil thermal, soil hydrology and carbon and nitrogen dynamics simulated in the TEM provided carbon substrate spatially and temporally for estimating the soil CO dynamics. Overall, although a few previous studies have examined the long-term impacts of climate, land use and nitrogen deposition on the CO dynamics (Chan and Steudler, 2006, Pihlatie et al., 2016), the global prediction of the soil CO dynamics still has a large uncertainty.

4.2 Major controls on soil CO dynamics

The sensitivity tests indicate that the consumption is normally much larger than the CO production, so that the former will determine the dynamics of the net flux (Table 5). The model being sensitive to air temperature explains the small increasing trends after the 1960s, the significant increasing trend in the 21st century and the large sinks over tropical areas (Table 5, Fig. 9). SOC did not directly influence the CO consumption. For instance, increasing SOC led to an increase in the soil CO substrate, implying that more CO in soils can be consumed. An extra 3 Tg CO yr⁻¹ was taken up from the atmosphere to the soil in the sensitivity test when SOC increased by 5 % (Table 5), which will be discussed in detail in Sect. 4.3. CO surface concentrations will only influence the uptake rate and soil CO substrate concentrations, thus influencing the soil CO consumption rate.

Table 6 Correlation coefficients between forcing variables (precipitation (Prec), air temperature (T_air), soil organic carbon (SOC), soil temperature (T_soil), soil moisture (M_soil) and atmospheric CO (CO air)) and absolute values of consumption, production and net flux for different regions and the globe.

The annual CO consumption and net flux have a similar correlation coefficient with forcing variables and both are significantly correlated with air temperature, soil temperature SOC and atmospheric CO concentration (R>0.91 globally, Table 6). Increasing temperature will increase microbial activities, while more SOC will increase soil CO substrate level. The annual CO consumption and net flux have low correlations with annual precipitation and soil moisture, especially at 45^∘ N–45^∘ S (R<0.54, Table 6). The annual CO production is strongly correlated with annual mean SOC, air temperature and soil temperature (R<0.91), while it is less correlated with precipitation, soil moisture and atmospheric CO concentration. Meanwhile, the monthly CO consumption, production and net flux are well correlated with air temperature, soil temperature, precipitation, and soil moisture (R>0.69 globally, Table 6). The soil moisture is significantly influenced by temperature at a monthly time step since the increasing temperature would induce higher evapotranspiration. The monthly CO consumption, production and net flux have low correlations with SOC because it will not change greatly within a month.

Figure 11 Global mean monthly time series of the MOPITT surface atmospheric CO concentration (ppbv) and soil CO consumption from model simulations E1 (Tg CO mon⁻¹).

The R between the annual soil CO consumption and atmospheric CO concentration is 0.91 at the global scale because the atmospheric CO concentration, air temperature, and soil temperature dominate the annual consumption rate. At a monthly scale, this R is −0.48 because the global atmospheric CO concentrations are high in winter and low in summer, while the simulated soil CO consumption shows an opposite monthly variation (Table 6, Fig. 11), suggesting that other factors such as precipitation, air temperature, and soil temperature are major controls for the monthly CO fluxes.

4.3 Model uncertainties and limitations

There are a number of limitations, contributing to our simulation uncertainties. First, due to lacking long-period observational data of the CO flux and associated environmental factors, the model parameterization can only be conducted for four ecosystem types, including boreal forest, temperate coniferous forest, temperate deciduous forest, and grassland. Tropical forest calibration is only conducted using a very limited amount of laboratory experiment data, but tropical areas are hotspots for CO soil–atmosphere exchanges. Besides, the amount of tropical forest SOC for the top 30 cm can be very large according to observations. The TEM may underestimate the top 30 cm SOC, which will underestimate the production rates, especially in tropical regions. Tropical regions typically have high temperature during the whole year, which may result in overestimation of the CO consumption using Eq. (2b). The large deviation of model simulations to observations in tropical savanna (which is a mosaic of tropical forest and grassland ecosystems) may be due to using outside air temperature to represent inside air temperature of transparent chamber observations (Varella et al., 2004), and uncertain tropical forest parameterization. Second, we used the conclusion from van Asperen et al. (2015) and only considered the thermal-degradation process for the CO production in this study. The photo-degradation process and biological formation process were not considered due to a lack of understanding of these processes. Third, the static CO surface concentration derived from the empirical function is lower than the MOPITT CO surface concentration, which will lead to underestimation of CO deposition velocity during 1901–2100. Fourth, from the sensitivity test (Table 5) we notice that an increase in SOC by 5 % resulted in a net flux increase from the atmosphere to the soil by 2.57 %. The SOC increase enhanced CO production (Eq. 3), CO concentrations (Eq. 1), and CO oxidation (Eq. 2). When the change in total oxidation is larger than the difference between the change in total production and the change in total soil CO concentration (Eq. 1), the estimate of the net flux change is negative (from the atmosphere to soil) using a mass balance approach (Sect. 2.2), which leads to a 2.57 % increase in the net flux in our SOC sensitive test. This is due to the fact that CO production (Eq. 3) is calculated independently from oxidation (Eq. 2). This will not influence our general findings since SOC varies only slightly during our simulation periods, with a 3 % increase from 1900 to 2013 (Fig. 4d) and up to a 4 % increase from 2014 to 2100 (Fig. 6g). This model artifact that is apparent in the SOC sensitivity test can be alleviated using a very fine time step (e.g., 1 s), because in this case CO concentrations change only slightly within the short time. Therefore, when a short time step is used, the net flux roughly equals the difference between production and oxidation. If the change in production is bigger than the change in oxidation, the change in net flux will be positive, leading to a decrease in deposition to the soil. The downside is that running the model at a time step of 1 s will require a significantly large amount of computing time. Fifth, our model structure still has a large potential to improve. In this study we divided the top 30 cm soil into 30 layers (layer thickness dz = 1 cm), but a finer division will increase the accuracy (Fig. 12). We chose 1 cm thickness because if thicker than 1 cm, the model vertical CO concentration profile will deviate from reality and the diffusion process will be influenced significantly. If thinner than 1 cm, it will need much more computing time but will not have much improvement compared to the thickness set to 1 cm (Fig. 12a–e). We notice that the 30-layer division represents the soil CO concentration profile well, not only in the days of soil CO net uptake, but also in the days of CO net emission (Fig. 12c, f). Sixth, the Michaelis–Menten function (Eq. 2a) is used in this model and we notice that k_CO is normally much larger than C(t,i) in those days of net soil uptake (over 10 times larger, Fig. 12). However, we cannot simplify Eq. (2b) to $f_{\mathrm{1}}\left(C\left(t,i\right)\right)=\frac{C\left(t,i\right)}{k_{\mathrm{CO}}}$, because the CO concentrations in the soil can be larger than in the atmosphere in the days of net emissions and C(t,i) may be close to k_CO, and then the simplified equation may lead to overestimation of CO oxidation (Fig. 12f). Finally, although we focused on natural ecosystems in this study, the land-use change, agriculture activity, and nitrogen deposition also affect the soil CO consumption and production (King and Crosby, 2002; Chan and Steudler, 2006). For instance, the soil CO consumption in agriculture ecosystems is from 0 to 9 mg CO m⁻² day⁻¹ in Brazil (King and Hungria, 2002). In this study, we used grassland or forest ecosystems to represent agriculture areas in the CODM module. Our future study will include these processes and factors.

Figure 12 Daily mean vertical soil CO concentration profiles of the top 30 cm. In the soil (depth < 0 cm), black diamonds represent the soil CO concentration (mg CO m⁻³). Above the surface (depth > = 0 cm), black diamonds represent the atmospheric CO concentration. (a), (b), (c), (d) and (e) are the results from the same day when soils are a net sink of CO, but using different layer thicknesses (dz = 10, 2, 1, 0.1 and 0.01 cm, respectively); (f) is the result from the day when soils are a net source of CO, with dz = 1 cm.

5 Conclusions

We analyzed the magnitude, spatial pattern, and controlling factors of the atmosphere–soil CO exchanges at the global scale for the 20th and 21st centuries using a process-based biogeochemistry model. Major processes include the atmospheric CO diffusion from the atmosphere to the soil and inside the soil of terrestrial ecosystems, microbial oxidation removal of CO, and CO production through chemical reaction. We found that air temperature and soil temperature play a dominant role in determining the annual soil CO consumption and production, while precipitation, air temperature, and soil temperature are the major controls for the monthly consumption and production. The atmospheric CO concentrations are important for annual CO consumption. We estimated that the global annual CO consumption, production and net fluxes for 2000–2013 are from −197 to −180, 34 to 36 and −163 to −145 Tg CO yr⁻¹, respectively, when using the MOPITT CO surface concentration data. Tropical evergreen forest, savanna and deciduous forest areas are the largest sinks, accounting for 66 % of the total CO consumption, while the Northern Hemisphere consumes 59 % of the global total. During the 20th century, the estimated CO deposition velocity is 0.16–0.19 mm s⁻¹. The predicted CO deposition velocity will reach 0.20–0.30 mm s⁻¹ in the 2090s, primarily because of the increasing air temperature. The areas near the Equator, eastern Asia, Europe and the eastern US will become the hotspots of sinks because they have warm and moist soils. This study calls for long-period observations of CO flux for various ecosystem types and better projection of atmospheric CO surface concentrations from 1901 to 2100 to improve future estimates of global soil CO consumption. The effects of land-use change, agriculture activities, nitrogen deposition, photo-degradation and biological formation will also be considered to improve future quantification of soil CO fluxes.