Temperature Sensitivity of Soil Respiration Probed by Numerical Analysis of Field-Observed Data Sets

Ippei Iiyama

School of Agriculture, Utsunomiya University, Utsunomiya, Tochigi, Japan.

DOI: 10.4236/gep.2023.118005                

Abstract

Temperature sensitivity of soil respiration is essential to predict possible changes in terrestrial carbon budget on various scenarios about atmospheric and soil climates. Although it is often evaluated by using respiratory quotient “Q₁₀”, Q₁₀ values of soil respiration seem to vary depending on methods or scales of evaluation. Aiming at probing how Q₁₀ values of soil respiration are evaluated differently for a field, this study used a model of soil respiration rate, and numerically evaluated soil respiration rates along depth by fitting the model to depth distributions of CO₂ concentration measured in a field. And temperature sensitivity of soil respiration rate was evaluated by comparing the determined soil respiration rates with atmospheric and soil temperatures measured in the field. The results showed that the relation between surface CO₂ emission rates and atmospheric temperatures was represented by lower Q₁₀ values than that between soil respiration rates and soil temperatures, presumably because the top soil layers had acclimatized in more extent to the existing thermal regime than the underlying deeper layers. Thus, for evaluating effects of long-term rise in atmospheric temperature on soil respiration, it is necessary to precisely predict the long-term change in depth distribution of soil temperature as well as to quantify temperature sensitivity of soil respiration along depth. The evaluated sensitivity of surface CO₂ emission rate to atmospheric temperature showed hysteresis, implying the needs for more knowledge about temperature sensitivity of soil respiration evaluated in both warming and cooling processes for better understandings and predictions about terrestrial carbon cycling.

Keywords

Air-Filled Porosity, Inverse Analysis, Mass Balance, Potentially Maximum CO2 Production Rate, Soil Gas Diffusion, Water Content

1. Introduction

Soil respiration is the sum of respiratory activities in a soil. Since respiratory activities in a soil must involve soil microbes or plant roots taking in free oxygen (O₂), decomposing organic substances, and producing metabolic energy, water, and carbon dioxide (CO₂), the intensity of the total of them is one of direct measures of soil health and can be quantified by measuring the amount of CO₂ released from the soil in a certain period.

Emission of CO₂ from a soil body attracts attention as an essential quantity concerning about carbon budgets between lands and atmosphere. This is because soils can be primary pool of carbon in terrestrial ecosystems. Past studies had estimated the size of global soil organic carbon to be up to 1395 [PgC] (Post et al., 1982) or 1576 [PgC] (Eswaran et al., 1993) . These kinds of estimations have been oddly updated with larger values year by year. For instance, Batjes (1996) evaluated the amount of global soil carbon within 1-meter depth as 2157 - 2293 [PgC]. Kochy et al. (2015) estimated it as nearly 3000 [PgC]. At the same time, the global soil respiration has been thought to emit 87 to 95 [PgC yr⁻¹] (Hashimoto et al., 2015) , which was equivalent to nearly ten times more amount of carbon emission than fossil-fuel burning and cement manufacture of 36.1 ± 0.3 [PgCO₂ yr⁻¹] (Liu et al., 2023) . These values suggest that changes in flow and storage of soil carbon may give large effects on the global carbon cycling.

Soil respiration is well related to regimes of ambient temperature and is commonly more promoted under warmer conditions in an ordinary temperature range. And the sensitivity of soil respiration to temperature has often been analyzed with the concept of respiratory quotient “Q₁₀” which is defined as a respiration rate at T + 10 [K] divided by that at T [K]. However, Q₁₀ values of soil respiration are likely to vary depending on methods or scales of evaluation. Relatively small Q₁₀ values of around 2 or less have been found mainly in global- or regional-scale studies of CO₂ emission from land surfaces, portion of which includes 1.5 by Knorr & Heimann (1995) , 1 to 2 by Kaminski et al. (2002) , 1.37 by Ise & Moorcroft (2006) , 1.43 to 2.03 by Zhou et al. (2009) , 1.5 by Bond-Lamberty & Thompson (2010) , 1.4 by Mahecha et al. (2010) , 1.41 to 1.43 by Zhan et al. (2017) .

However, CO₂ emission from soils in site-specific or soil-sample scales seems to be characterized with relatively larger Q₁₀ values of more than 2. For example, Raich & Schlesinger (1992) compiled the Q₁₀ values from 15 literatures and estimated the median value of 2.4. Wang et al. (2010) also collected 185 data sets from 69 publications, and obtained the average of 2.67. Wang et al. (2018) measured soil respiration rates in a meadow steppe for two years, and their seasonal and annual Q₁₀ values ranged from 2 to 4. Kim et al. (2020) also measured soil respiration rates throughout a year in an urban and a well reserved forest, and reported Q₁₀ values of more than 3.

These apparent differences among the estimated Q₁₀ values has not yet been explained well, and may cause such a hypothesis that a Q₁₀ value of soil respiration is estimated higher when soil temperatures are related to respiration rates in the soil layer than when atmospheric temperatures are adopted to explain CO₂ emission rates from the surface of the soil layer, since the studies about soil respiration in larger spatial scales often monitor above-ground meteorology while the studies in soil-sample scales or local-site scales had often correlated respiration rates with temperatures of their soils.

One promising way to examine this hypothesis is to use a mathematical model that allows to evaluate both a surface CO₂ emission rate from a soil layer and a depth distribution of soil respiration rates in the soil layer for given sets of measured data including gaseous CO₂ concentration in soil and temperatures of both the soil layer and the atmosphere. Thus, this study proposed a model of soil respiration rate, and determined depth distribution of soil respiration rates by making the model’s outputs best reproduce the dynamics of soil gaseous CO₂ measured in a field. And, for examining the hypothesis mentioned above, the determined soil respiration rates were compared with both atmospheric and soil temperatures monitored in the field.

2. Materials and Methods

2.1. Field Data to Run the Model

The inputs for the model of soil respiration rate were the data sets observed by Iiyama (2023) in a meadow field. The field was underlaid with volcanic ash soil layers which consisted of a topsoil layer with the thickness of about 0.25 [m], followed by a transitional layer of about 0.3 [m] in thick, underlain by a subsoil layer. Two study sites were set to represent a tree-standing area and a harvesting area in the field, and named the sites A and B, respectively. The range of rooting depth in the site A reached at least 1 [m] in depth, while that in the site B was found mainly within 0.25 [m] from the soil surface.

The data sets consisted of the depth profiles of state variables and physical properties of the soil layers in the field. The group of the state variables was comprised of CO₂ concentration of soil gas c [kg·m⁻³], soil temperature T_s [K], and volumetric water content θ [m³·m⁻³], which were measured at several depths each as shown in Figure 1. The soil physical properties included soil bulk density ρ_d [kg·m⁻³] and soil particle density ρ_s [kg·m⁻³] (Table 1). By using these two

Table 1. The soil bulk densities, particle densities, and saturated volumetric water contents of the two sites (Data from Iiyama (2023) ). Each value is an average of triplicate measurements with a standard deviation in the parentheses. The soil textures of all the soil layers were classified as clay loam soils on the basis of the soil-texture classification defined by the International Union of Soil Science (IUSS).

Figure 1. The time-series data sets of CO₂ concentration in soil air (i-A and i-B), volumetric water content (ii-A and ii-B), and atmospheric and soil temperatures (iii-A and iii-B) (Data from Iiyama (2023) ). Either of “A” and “B” in the sub-indices of the sub-graph labels denotes the site name where the data set on the chart was measured in the field. The numerals in a sub-graph indicate the depths of measurement in centimeter. The values of temperature are plotted on a Celsius scale, instead of a Kelvin scale.

physical quantities, saturated volumetric water content θ_s [m³·m⁻³] was evaluated as 1 − ρ_d/ρ_s. The details about the measurements of these state variables with the physical properties were explained in Iiyama (2023) . In parallel with the state variables, atmospheric temperature T_a [K] at each site was also measured by a temperature logger (HOBO U23 Pro v2 external temperature/relative humidity data logger (U23-002); Onset Computer Corp; Bourne, MA, USA), which was contained in a solar radiation shield (RS1; Onset Computer Corp; Bourne, MA, USA).

2.2. Modeling Soil Respiration Rate

When the respiratory activity in a soil is represented by an amount of CO₂ produced in the soil, some optimum condition can be supposed, in which a unit mass of soil emits CO₂ in an optimum rate s_o [kg·kg⁻¹·s⁻¹]. When the soil has a bulk density of ρ_d [kg·m⁻³], the quantity s_o can be converted into the optimum rate of CO₂ production per unit volume of the soil ρ_ds_o [kg·m⁻³·s⁻¹].

Since the respiratory activity requires sure supply of oxygen, it should be restricted when soil moisture content is too high to allow soil air to be smoothly exchanged. At the same time, the respiratory activity weakens under an extremely dry condition because of the lowering of metabolic activities of soil microbes or plant roots with the lack of soil moisture causing the inhibition of flows of substrate supply to the activities. Thus, for modeling the soil respiration rate under a field condition, the term ρ_ds_o should be multiplied by some inhibiting factor ι [non-dim] which takes a value between 0 and 1, depending on water content θ [m³·m⁻³] or air-filled porosity a [m³·m⁻³] of the soil. Therefore, by using the expressions s_o, ρ_d, and ι, the soil respiration rate S [kg·m⁻³·s⁻¹] can be modeled with the following expression:

$S=\iota {\rho }_d{s}_o$ (1)

According to past studies for soil respiration in fields (Liu et al., 2002; Xu et al., 2004; Jassal et al., 2008; Zhang et al, 2010) and in laboratories (Guntinas et al, 2013; Zhang et al, 2015) , soil respiration is likely to be limited in the lowest and highest extremities of the entire moisture range. Thus, this study modeled ι as below:

$\iota =\frac{a}{a+{\kappa }_0}\frac{\theta }{\theta +{\kappa }_0}{\kappa }_1=\frac{{\theta }_s-\theta }{{\theta }_s-\theta +{\kappa }_0}\frac{\theta }{\theta +{\kappa }_0}{\kappa }_1$ (2)

where κ₀ [m³·m⁻³] is a constant that regulates the ways of increase and decrease in ι with the changes in θ. The other parameter κ₁ can be mathematically derived when it is set in order for making ι take 1 as its maximum as follows:

${\kappa }_1={\left(\frac{{\theta }_s}{{\theta }_s+2{\kappa }_0}\right)}^{-2}$ (3)

The optimum CO₂ production rate s_o in Equation (1) should vary along depth due to the differences in population densities of microbes or plant roots in a soil layer, meaning that s_o can be modeled as a function of the vertical location. As a first approximation, it was assumed that there is a depth at which s_o takes the maximum along depth. As an example of smooth curve that satisfies the above assumption, this study used such an expression as:

$s_o=\{\begin{array}{l}0\left(z_u\le z\le 0\text{or}z\le z_l\right)\\ \frac{s_{o\max }}{2}\left(1+\cos \left(\pi \frac{z-z_{\max }}{z_u-z_{\max }}\right)\right)\left(z_{\max }\le z\le z_u\right)\\ \frac{s_{o\max }}{2}\left(1+\cos \left(\pi -\pi \frac{z-z_l}{z_{\max }-z_l}\right)\right)\left(z_l\le z\le z_{\max }\right)\end{array}$ (4)

where z [m] indicates the upward-positive vertical location and takes 0 [m] at the level of a soil surface, z_u [m] and z_l [m] are the upper most and the lower most levels of the domain of soil respiration, s_{o max} [kg·kg⁻¹·s⁻¹] is the potentially maximum rate of CO₂ production in the domain $z_l\le z\le z_u$ , and z_max [m] is the vertical location where s_o takes s_{o max}.

The model described as Equations (1) to (4) means that the depth profile of soil respiration rate S is explained by the soil physical state that the depth distributions of ρ_d, θ, and θ_s represent and by such respiratory characteristics of the soil layer as s_{o max}, z_l, z_max, z_u, and κ₀. Among these expressions, s_{o max} is the only expression that can be related to soil temperature T_s [K]. Therefore, for probing the temperature sensitivity of S, s_{o max} and z_max were set as unknown variables to be identified for given sets of the other terms. And the identified series of s_{o max} was compared with a series of T_s obtained at corresponding times and depths.

The rest of the parameters were evaluated by using the measured data sets. ρ_d and θ_s in Equations (1) to (3) were given by citing the values tabulated on Table 1. The time-series data of θ on Figure 1 was used to assign values for both θ and a in Equation (2) and (3). z_u = 0 [m] and z_l = −1 [m] were assigned with the assumption that respiration activity may be found at any depth in the soil layer of analysis. κ₀ = 0.1 was assigned so that ι decreases sharply in the top-most and the bottom-most 10- to 20-percent of the whole domain of θ. The way to identify s_{o max} and z_max with these settings will be explained in the following sections.

2.3. Mass Balance about CO₂ in a Soil Layer

To set a problem for identifying s_{o max} and z_max in Equation (4) is to set up a series of equations that involves S. Thus, the mass balance about CO₂ in a soil layer was formulated for the development of the equations. As the first step of the application, a soil layer to be analyzed was discretized with a series of calculation nodes arranged along depth. The nodes were numbered from 0 to n to point any of vertical locations in the soil layer from the bottom z₀ [m] to the surface z_n [m]. Then, gaseous CO₂ concentrations at z = z_i ( $0\le i\le n$ ) can be denoted as c_i [kg·m⁻³] on the node i. When a certain c_i is needed to be compared with a measured CO₂ concentration, the measured data of CO₂ concentration in Figure 1 were spatially interpolated and used.

A segment was defined for every node so that c_i represents the value of CO₂ concentration for the segment that includes the node i. Any of the segment boundaries was placed at the center of two adjacent nodes so that a segment around z = z_i spans between (z_i + z_i−1)/2 and (z_i+1 + z_i)/2 for $0\le i\le n-1$ . For i = 0 and i = n, the segments cover the domains $z_0\le z\le \left(z_{\text{1}}+z_0\right)/2$ and $\left(z_n+z_n{}_{-1}\right)/2\le z\le z_n$ , respectively. In this study, n = 20, z₀ = −1 [m], and z_n = z₂₀ = 0 [m] were assigned, and all of z_i for $\text{1}\le i\le n-\text{1}$ were arranged every 0.05 [m] along depth, for the soil layers of analysis was 1 [m] in thick and most of the measured data sets were obtained with the spatial resolution of more than 0.1 [m] as shown in Figure 1 and Table 1.

Then, the flow of gaseous CO₂ in the soil layer was formulated. Since gas flow in soils owes almost solely on gaseous diffusion driven by difference in concentration of the gas species along depth, a value of CO₂ gas flux q_i [kg·m⁻²·s⁻¹] can be evaluated on a segment boundary between z = z_i and z = z_i+1 such that:

$q_i=-D_{si}\frac{c_{i+1}-c_i}{z_{i+1}-z_i}$ (5)

where D_{s i} [m²·s⁻¹] is the soil gas diffusion coefficient for CO₂ assigned to the domain covering between the centers of two adjacent segments i and i + 1. The formulation of D_s will be described in Section 2.4.

By using the terms defined above, the mass balance equations of CO₂ for the segments i ( $\text{1}\le i\le n-\text{1}$ ) and two boundary conditions can be expressed as follows:

$\{\begin{array}{ll}\frac{M_{i,j+1}-M_{i,j}}{t_{j+1}-t_j}\frac{z_{i+1}-z_i}{2}=q_b-q_i+S_i\frac{z_{i+1}-z_i}{2}\hfill & \left(i=0\right)\hfill \\ \frac{M_{i,j+1}-M_{i,j}}{t_{j+1}-t_j}\frac{z_{i+1}-z_{i-1}}{2}=q_{i-1}-q_i+S_i\frac{z_{i+1}-z_{i-1}}{2}\hfill & \left(1\le i\le n-1\right)\hfill \\ c_{i,j+1}=c_a\hfill & \left(i=n\right)\hfill \end{array}$ (6)

where M_i,j [kg·m⁻³] is the mass of CO₂ per unit volume of soil stored in the segment i at a time t_j [s], being called the storage term in this study. The subscript j denotes a time-step number (j ≥ 0) so that the left sides of the first and second sub-equations in Equation (6) denote the rates of increase in the mass of CO₂ stored in a segment i for a time interval Δt = t_j+1 − t_j [s]. In this study, 1800[s] was assigned to Δt as an increment sufficiently small compared with the measurement intervals for the state variables c, T_s, and θ in the field (Figure 1). The formulation about M_{i, j} will be explained in section 2.5. S_i [kg·m⁻³·s⁻¹] is the mass of CO₂ produced per unit volume of soil per unit time in the segment i, and was modeled as Equation (1). q_b [kg·m⁻²·s⁻¹] is the CO₂ gas flux flowing into the segment i = 0 from the bottom of the segment, and, c_a [kg·m⁻³] is the atmospheric CO₂ concentration. Both q_b and c_a were given as the lower and the upper boundary conditions of the spatial domain of analysis. In this study, 0 [kg·m⁻²·s⁻¹] was assigned to q_b as the zero-flux lower boundary condition for imitating the CO₂ concentration gradient at around z = −100 [cm] in the field having been almost 0 during the monitoring period (Figure 1). This study also assigned a value equivalent to 0.04 [%] to c_a as the atmospheric upper boundary condition. All of the expressions q_b, q_i, q_i−1, S_i, and c_a in the right sides of the three sub-equations in Equation (6) were defined as their time-averaged values between t = t_j and t = t_j+1 so that the Crank-Nicholson scheme was applied to the process of solving Equation (6).

2.4. Soil Gas Diffusion Coefficient

Soil gas diffusion coefficient D_{s i} [m²·s⁻¹] in Equation (5) was evaluated as follows:

$D_{si}={\xi }_i{D}_{ai}$ (7)

where ξ_i [non-dim] is the pore tortuosity factor while D_{a i} [m²·s⁻¹] is the diffusion coefficient of gaseous CO₂ in the atmosphere. These two terms were both assigned to each of the domains of $z_i\le z\le z_i{}_{+\text{1}}$ ( $0\le i\le n-\text{1}$ ).

The tortuosity factor ξ_i in Equation (7) is a function of air-filled porosity of a soil, and was expressed by the following model (Iiyama et al., 2012) :

${\xi }_i=\{\begin{array}{ll}0\hfill & \left(\bar{a}<a_0\right)\hfill \\ \lambda {\left(\bar{a}-a_0\right)}^{m_1}/{\bar{\theta }}_s^{m_2}\hfill & \left(\bar{a}\ge a_0\right)\hfill \end{array}$ (8)

where $\bar{a}$ [m³·m⁻³], $\bar{\theta }$ [m³·m⁻³], and ${\bar{\theta }}_s$ [m³·m⁻³] are the spatially-averaged air-filled porosity, volumetric water content, and saturated volumetric water content of the soil which are evaluated by citing θ and θ_s at z = z_i and z_i+1, a₀ [m³·m⁻³] is the minimum air-filled porosity for gaseous diffusion practically taking place in the soil, λ [non-dim], m₁ [non-dim], and m₂ [non-dim] are the constants that vary with type of soil. The values of the parameters applied in this study are tabulated in Table 2, which were obtained by fitting Equation (8) to the data sets of soil gas diffusivity for Andosol layers in the study field measured by Iiyama (2023) .

The diffusion coefficient of gaseous CO₂ in the atmosphere, D_{a i} in Equation (7), was expressed as a function of pressure and temperature around the place of interest as follows (Freijer & Leffelaar, 1996; Osozawa & Hasegawa, 1995; Makita, 1988) :

$D_{ai}=D_{astd}{\left(\frac{{\bar{T}}_s}{T_{sstd}}\right)}^{\nu }\frac{p_{std}}{\bar{p}}$ (9)

where D_{a std} = 1.38 × 10⁻⁵ [m²·s⁻¹], T_{s std} = 273.15 [K], ν = 1.75, and p_std = 101.325 [kPa]. ${\bar{T}}_s$ [K] is the soil temperature in the domain of $z_i\le z\le z_i{}_{+\text{1}}$ which is evaluated by citing T_s at z = z_i and z_i+1. For the pressure of the soil air, p_std was assigned to $\bar{p}$ with the assumption that most of the air-filled pores in the entire soil layer are opened to the atmosphere.

2.5. CO₂ Storage Term

The storage term M_{i, j} [kg·m⁻³] in Equation (6) was expressed as the sum of the amount of CO₂ contained in both the liquid and gaseous phases in a segment i,

Table 2. The parameters of the tortuosity ξ in Equation (8) for the soil layers in the study field (Data from Iiyama (2023) ).

since the solubility of CO₂ is relatively high among gas species in soil air. When the dissolved fraction of CO₂ is considered, the storage term can be expressed as:

$M_{i,j}=\left(a_{i,j}+{\theta }_{i,j}{K}_H\right)c_{i,j}=\left({\theta }_{si}-{\theta }_{i,j}+{\theta }_{i,j}{K}_H\right)c_{i,j}$ (10)

where a_{i, j} and θ_{i, j} are air-filled porosity and volumetric water content of soil in the segment i at a time t = t_j. The Henry’s coefficient for CO₂, K_H, was introduced with the assumption that gas-liquid equilibrium can be achieved instantly. For K_H, the following polynomial function was used for considering its temperature dependence:

$K_H=k_0+k_1{T}_{si,j}+k_2{T}_{si,j}^2+k_3{T}_{si,j}^3$ (11)

where T_{s i,j} is soil temperature in the segment i at a time t = t_j. The parameters k₀ = 217.8098, k₁ = −1.990097, k₂ = 6.096563 × 10⁻³, and k₃ = −6.250000 × 10⁻⁶ were obtained by fitting the third-order polynomial function to the reference values tabulated on National Astronomical Observatory of Japan (2020) .

2.6. Evaluation of Soil Respiration Rates

The numerical scheme made of the series of equations from (1) to (11) can be used for obtaining a time-series of CO₂ concentrations c_i along depth for a given set of state variables, physical properties, and parameters including s_{o max} and z_max. This study used this numerical scheme inversely to identify an optimum pair of s_{o max} and z_max by fitting the numerically-determined c_i to the values of CO₂ concentration obtained in the field in the following steps (Figure 2):

Step 1 Select a date of measurement of soil gaseous CO₂ concentrations in the field from the time domain plotted in Figure 1 and name it t₀.

Step 2 Select the date of measurement next to t₀, and name it t₁. Since the interval of CO₂ concentration measurements was 7 days, t₁ - t₀ was equivalent to 7 days.

Step 3 Assign the measured values of soil gaseous CO₂ concentrations at the two consecutive measurements to the sets of variables c_meas(z, t₀) and c_meas(z, t₁).

Step 4 Let c_meas(z, t₀) be the initial condition for the numerical scheme.

Step 5 Set up the input data sets for the numerical scheme. The inputs are comprised of the state variables other than CO₂ concentrations, the physical properties, and the parameters except for s_{o max} and z_max. The values of θ and T_s for any time and any location were evaluated by linearly interpolating the measured values obtained between t₀ and t₁ at two neighboring measurement depths.

Step 6 Define the domains and increments about s_{o max} and z_max for finding out an optimum pair $s_{o\max }^{\ast }$ and $z_{\max }^{\ast }$ that minimizes the error between a numerically-determined c_i, namely c_calc(z, t₁), to the measured values c_meas(z, t₁).

Step 7 Prepare a tentative pair of s_{o max} and z_max from the domains defined in Step 6.

Step 8 Run the numerical scheme and obtain a time-series of c_i as the output for the time domain of analysis being $0\le t\left[\text{s}\right]\le t_{\text{max}}$ , where t_max = 7 × 86400 [s].

Step 9 Assign the calculated values of soil gaseous CO₂ concentrations at t = t_max to the set of variables c_calc(z, t₁).

Figure 2. The flow chart to outline the process for identifying an optimum pair of s_{o max} and z_max by fitting a numerically determined CO₂ concentration profile to a measured CO₂ concentration profile.

Step 10 Evaluate the error between c_meas(z, t₁) and c_calc(z, t₁) with the concept of least squares.

Step 11 If the error is smaller than the minimum of ever-recorded errors, then update a best-ever pair $s_{o\max }^{\ast }$ and $z_{\max }^{\ast }$ with the currently-tried s_{o max} and z_max, and the minimum-ever error is also updated with the error just evaluated in Step 10.

Step 12 If the domain defined in step 6 has not yet entirely been scanned, return to step 7. If it has been done, exit the iteration of the series of steps from Step 7 to this Step 12 and obtain the optimum pair of s_{o max} and z_max as $s_{o\max }^{\ast }$ and $z_{\max }^{\ast }$ .

By applying the series of these steps to every measurement interval of soil gaseous CO₂ concentration in the field, a time-series of s_{o max} and z_max were determined.

After the determination of a time-series of s_{o max} and z_max, time-series of CO₂ emission rates from the soil surfaces q_sur [kg·m⁻²·s⁻¹] were also evaluated for all the measurement intervals by applying Equation (5) to the top-5-cm layer and using c_calc(z, t) for the identified pairs of s_{o max} and z_max. Thus, q_sur is expressed as follows:

$q_{sur}=\frac{1}{J}\underset{j=1}{\overset{J}{{\displaystyle \sum }}}\left(-D_s\frac{c_{\text{calc}}\left(z_n,t_j\right)-c_{\text{calc}}\left(z_{n-1},t_j\right)}{z_n-z_{n-1}}\right)$ (12)

where J is the maximum of the time-step number j such that t_J − t₀ is equivalent to a measurement interval of CO₂ concentration profile, being 7 days in this study. The values of c_calc(z_n, t_j) and c_calc(z_n−1, t_j) are the CO₂ concentrations calculated for the time t_j and the vertical locations of z = 0 [m] and z = −0.05 [m], respectively. D_s was evaluated in the same manner as explained in section 2.4.

3. Results

Figure 3(i) and Figure 3(ii) show the time-series of s_{o max} and z_max in Equation (4), identified by the numerical scheme explained in Section 2.6. The scale of s_{o max} is in “ngCO₂ g⁻¹ s⁻¹”. The values of s_{o max} increased from May to August, decreased by December and, then, remained almost constant in their lowest levels till the next spring. The maximum values of s_{o max} were recorded in the end of July with 2 [ngCO₂ g⁻¹ s⁻¹] in the site A and 0.88 [ngCO₂ g⁻¹ s⁻¹] in the site B, respectively, while the lowest values in the winter term were 0.22 and 0.08 [ngCO₂ g⁻¹ s⁻¹] in the sites A and B, respectively. The values for the site A were two times or more of those for the site B in any time, suggesting that respiratory activity under the tree-standing area was basically stronger than that under the harvesting area.

The seasonal trends of z_max indicated that the most active depth ranges of soil respiration stayed in the top soil layers from May to September, while they went down to the subsoil layers during the winter term. These evaluations can largely be explained by the seasonal change in depth distribution of soil temperature,

Figure 3. The time-series of (i) the potentially maximum rate of CO₂ production along depth s_{o max}, (ii) the vertical location z_max where s_{o max} emerges, and (iii) surface CO₂ emission rate q_sur.

implying that the depth range of soil respiration would follow the seasonal transition of the warmest depth in a soil profile. However, there was also two-month difference in the commencement of deepening in z_max between the two sites. This difference was presumably because the respiratory activities of the larch tree roots in the site A have more significant tolerance to the lowering of soil temperature than those of the grasses and forbs in the site B.

Figure 3(iii) shows the seasonal behaviors of surface CO₂ emission rate q_sur evaluated by using Equation (12). The emission rates ranged between 17 and 125 [MgCO₂ ha⁻¹ yr⁻¹] in the site A, while between 5 and 61 [MgCO₂ ha⁻¹ yr⁻¹] in the site B, being almost equivalent to the values reported in past studies about carbon emission from grassland fields (Risk et al., 2002; Fierer et al., 2005; Lee et al., 2007) .

Some local minimum values were found in the midst of June and in the end of September while z_max stayed in the top soil layers. Each of them was presumably due to sudden momentary lowering of atmospheric temperature (Figure 1(iii)). In these situations, soil temperature in most part of the top soil layers also fell clearly, accompanying with the decrease in atmospheric temperature. These facts indicate the quickness and sureness of the response of surface CO₂ emission to the change in atmospheric temperature, which can be found particularly when the depth range of most active respiration resides near the soil surface.

On the other hand, in August, momentary reduction in surface CO₂ emission coincided with z_max deepening into the subsoil layers. In the midst of August, there was a dry period lasting more than two weeks, and soil water contents in the top soil layers dropped to its lowest level of the year (Figure 1(ii)). This dried condition was a possible cause of the momentary reduction in soil respiratory activity, and it was implied that soil moisture condition can also change the most active depth range of soil respiration while soil thermal condition can primarily affect the intensity of soil respiration at each depth.

4. Discussion

Figure 4 shows the potential maximum rates of CO₂ production s_{o max} [ngCO₂ g⁻¹ s⁻¹] plotted against soil temperature in Celsius scale. Each of the two sub-graphs includes the data set obtained from each site. And the data set in each sub-graph is divided into three sub-groups, each of which represents a certain depth range among the top, transitional, and subsoil layers as classified in Table 1 and Table 2.

In both sites, the values for the top soil layers were likely to be higher than those for the other deeper layers for a given temperature, reflecting the denser populations of plant roots and soil microbes in the shallower depth ranges.

The rates of increase in s_{o max} tended to increase with rising temperature regardless of depth range, suggesting that the relationship between soil temperature T_s [K] and s_{o max} [ngCO₂ g⁻¹ s⁻¹] can be described as follows:

Figure 4. The temperature dependence of the potentially maximum rate of CO₂ production along depth s_{o max} [ngCO₂ g⁻¹ s⁻¹]. The sub-graphs (i) and (ii) include the values for the sites A and B, respectively. On each sub-graph, soil temperature T_s is indicated on a Celsius scale, instead of a Kelvin scale. The numerals in the legend of each sub-graph denote the depth ranges which the sub-groups of the data set represent each, described in meter.

$s_{o\max }={\alpha }_s\exp \left({\lambda }_s\left(T_s-273.15\right)\right)$ (13)

Thus, Equation (13) was fitted to each of the data sets on Figure 4. The results were tabulated on Table 3 with the two parameters α_s [ngCO₂ g⁻¹ s⁻¹] and λ_s [K⁻¹] identified through the curve-fitting processes. The values of Q₁₀ on Table 3 were calculated by using such an expression as Q₁₀ = exp(10λ_s) in accordance with its definition.

Equation (13) reproduced well the obtained T_s − s_{o max} relations. Since the soil respiration model in this study had no a-priori assumption about temperature sensitivity, it was inferred that the relation between soil temperature and potentially maximum CO₂ production rate in the study field can be commonly characterized by the exponential-type function.

The temperature sensitivity of soil respiration rate differed among the depth ranges. In both sites, the top soil layers showed the lowest Q₁₀ values in the entire soil layers. This means that respiratory activities in the deeper soil layers can be more easily stimulated by the change in soil thermal regime, while the absolute amounts of CO₂ production were larger in the top soil layers than in the deeper layers.

In summation, the relatively simple model of soil respiration combined with in-situ monitoring of CO₂ concentration profile served as an effective measure for both identifying the depth range of the most active layer of soil respiration and characterizing the temperature sensitivity of soil respiration in each soil layer.

Figure 5 shows the temperature dependence of the surface CO₂ emission rate q_sur. The horizontal axis denotes the atmospheric temperature in Celsius scale. The two sub-graphs are served for contrasting the two data sets obtained from the sites A and B. The data set in each sub-graph is divided into six sub-groups to describe the temporal change in the T_a − q_sur relations. The sub-groups are labelled with month-numbers from May and June with 5 and 6 to March with 3.

The T_a − q_sur relations in both sites were also characterized as follows:

$q_{sur}={\alpha }_q\exp \left({\lambda }_q\left(T_a-273.15\right)\right)$ (14)

Table 3. The parameters about the temperature dependence of soil respiration rate, evaluated by fitting Equation (13) to each of the data sets of T_s − s_{o max} relations on Figure 4. The fitting process is basically a linear regression in which Equation (13) is transformed into such a linear expression as “log_es_{o max} = log_eα_s + λ_s(T_s − 273.15)”. R² is the coefficient of determination of the linear regression. Q₁₀ is calculated in accordance with its definition “exp(10λ_s)”.

Figure 5. The temperature dependence of the surface CO₂ emission rate q_sur [MgCO₂ ha⁻¹ yr⁻¹]. The sub-graphs (i) and (ii) describe the data plots for the sites A and B, respectively. On each sub-graph, atmospheric temperature T_a is indicated on a Celsius scale, instead of a Kelvin scale. The numerals in the legend of each sub-graph denote the month numbers for which the corresponding data set was obtained.

And fitting Equation (14) to the values on Figure 5 gave the values of the two parameters α_q [MgCO₂ ha⁻¹ yr⁻¹] and λ_q [K⁻¹] as shown in Table 4.

The obtained Q₁₀ values for T_a − q_sur relations were 1.968 and 2.153 for the sites A and B, respectively. The comparisons of these values with the Q₁₀ values for the transitional and the subsoil layers listed on Table 3 suggest that the sensitivity of surface CO₂ emission rate to atmospheric temperature can be represented by lower Q₁₀ values compared with the sensitivity of soil respiration rate to temperature at the depth of CO₂ production, supporting the hypothesis of this study.

On the other hand, these values were similar in size to the Q₁₀ values for the top soil layers of the two sites, implying that the respiratory activities in the top soil layers were main source of the surface CO₂ emission during the field measurement period, and had acclimatized in more extent to the existing thermal regime than the other deeper layers so that they were relatively insensitive to change in soil temperature. These implications mean that the seasonal variation of surface CO₂ emission rate from a soil layer can be approximately predicted by evaluating the temperature sensitivity of the surface CO₂ emission rate itself or of the CO₂ production rate near the surface of the soil layer. At the same time, however, it is also clarified that for evaluating possible effects on soil respiration of long-term rise in atmospheric temperature through which a yearly-mean soil temperature can change, it is necessary to precisely predict the long-term change in depth distribution of soil temperature as well as to quantify temperature sensitivity of soil respiration at each depth.

The T_a − q_sur relations showed seasonal hysteretic behaviors (Figure 5). The size of hysteretic loop on the T_a − q_sur relations was larger for the data sets about the site A than about the site B. In the warming season, q_sur is likely to take smaller

Table 4. The parameters about the temperature dependence of surface CO₂ emission rate, evaluated by fitting Equation (14) to each of the data sets of q_sur on Figure 5. The fitting process is basically a linear regression in which Equation (14) is transformed into such a linear expression as “log_eq_sur = log_eα_q + λ_q(T_a − 273.15)”. R² is the coefficient of determination of the linear regression. Q₁₀ is calculated in accordance with its definition “exp(10λ_q)”.

values for a given temperature than in the cooling season. On Figure 5(i), the temporal averages of q_sur in the domain of atmospheric temperature from 15 to 25 [˚C] were 69 [MgCO₂ ha⁻¹ yr⁻¹] between May and June and 102 [MgCO₂ ha⁻¹ yr⁻¹] between September and October. Similarly, the values of q_sur for the November-December period were also clearly higher than those for the January-March period in the domain of T_a between −5 and 10 [˚C]. The T_a − q_sur relations for the site B also showed hysteresis in almost the same pattern as the results for the site A. But the size of the hysteretic loop for the site B was very small. These tendencies imply that the hysteretic behavior in T_a − q_sur relations is likely to become more significant for a soil layer with stronger respiratory activities, and that the temperature sensitivity of soil respiration may be over- or under-estimated if it is assumed to be a single-valued function. Thus, it is recommended to gather more knowledge about temperature sensitivity of soil respiration obtained in both warming and cooling processes for better understandings and predictions about terrestrial carbon cycling.

5. Conclusion

This study proposed a model of soil respiration rate for determining depth distribution of soil respiration rates based on the dynamics of soil gaseous CO₂ measured in a field. Then, temperature sensitivity of soil respiration rate was evaluated by using both atmospheric and soil temperatures in the field, with such a hypothesis that a Q₁₀ value of soil respiration is estimated higher when respiration rates in a soil layer are related to soil temperatures than when surface CO₂ emission rates from the soil layer are explained by atmospheric temperature.

The numerical scheme in this study clarified the seasonal behavior of the most active depth range of soil respiration, and implied that the depth range of soil respiration follows the seasonal transition of the warmest depth in a soil profile, though it seems also necessary to consider some plant physiological features like the difference in temperature-tolerance among plant species.

The numerical scheme in this study reproduced annual trends of surface CO₂ emission rate with acceptable values, suggesting that a sudden momentary change in atmospheric temperature directly affects surface CO₂ emission rate without changing the most active depth range of soil respiration, while an extremely dry condition can change, though tentatively, most active depth range of soil respiration itself.

The relation between soil temperature and potentially maximum CO₂ production rate obtained in this study indicated that the temperature sensitivity of soil respiration rate can be reproduced well by the exponential-type function. According to the obtained relation, the top soil layers had the lowest Q₁₀ values in the entire soil layers, meaning that respiratory activities in deeper soil layers can be more sensitive to the change in soil thermal regime, even though the absolute amounts of CO₂ production are larger in the top soil layers so far.

The sensitivity of surface CO₂ emission rate to atmospheric temperature was characterized by lower Q₁₀ values compared with that of soil respiration rate to soil temperature, supporting the hypothesis of this study. Therefore, it was inferred that for evaluating possible effects on soil respiration of long-term rise in atmospheric temperature through which a yearly-mean soil temperature can change, it is necessary to precisely predict the long-term change in depth distribution of soil temperature as well as to quantify temperature sensitivity of soil respiration at each depth.

However, it was also noticed that the main source of surface CO₂ emission in a field is often in the top soil layer of the field, where respiratory activities may have acclimatized to the existing thermal regime. Thus, seasonal variation of surface CO₂ emission rate from a soil layer can be approximately predicted by evaluating the temperature sensitivity of the surface CO₂ emission rate itself or of the CO₂ production rate near the surface of the soil layer.

The evaluated sensitivity of surface CO₂ emission rate to atmospheric temperature showed hysteresis. And the tendency was more significant for a soil layer with larger potentially maximum CO₂ production rate. This means that when the temperature sensitivity of soil respiration is modeled by a single-valued function, the model may over- or under-estimate actual soil respiration rates. Therefore, further studies are required to obtain more knowledge about temperature sensitivity of soil respiration in both warming and cooling processes for better understandings and predictions about terrestrial carbon cycling.

Acknowledgements

The author thanks Mr. T. Shiozawa, Mr. E. Saito and Mr. N. Yamaguchi in the Utsunomiya University Farm for their supports in the managements of the study field. The author also thanks Mr. Suto in the Graduate School of Agriculture, Utsunomiya University, and Ms. M. Anzai, Mr. K. Ogasawara, Mr. Y. Sakanishi, and Mr. Y. Usui in the School of Agriculture, Utsunomiya University, for their measurement works about the soils of the study field.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Batjes, N. H. (1996). Total Carbon and Nitrogen in the Soils of the World. European Journal of Soil Science, 47, 151-163. https://doi.org/10.1111/j.1365-2389.1996.tb01386.x
[2]	Bond-Lamberty, B., & Thomson, A. M. (2010). Temperature-Associated Increases in the Global Soil Respiration Record. Nature, 464, 579-582. https://doi.org/10.1038/nature08930
[3]	Eswaran, H., van den Berg, E., & Reich, P. (1993). Organic Carbon in Soils of the World. Soil Science Society of America Journal, 57, 192-194. https://doi.org/10.2136/sssaj1993.03615995005700010034x
[4]	Fierer, N., Chadwick, O. A., & Trumbore, S. E. (2005). Production of CO2 in Soil Profiles of a California Annual Grassland. Ecosystems, 8, 412-429. https://doi.org/10.1007/s10021-003-0151-y
[5]	Freijer, J. I., & Leffelaar, P. A. (1996). Adapted Fick’s Law Applied to Soil Respiration. Water Resources Research, 32, 791-800. https://doi.org/10.1029/95WR03820
[6]	Guntinas, M. E., Gil-Sotres, F., Leiros, M. C., & Trasar-Cepeda, C. (2013). Sensitivity of Soil Respiration to Moisture and Temperature. Journal of Soil Science and Plant Nutrition, 13, 445-461. https://doi.org/10.4067/S0718-95162013005000035
[7]	Hashimoto, S., Carvalhais, N., Ito, A., Migliavacca, M., Nishina, K., & Reichstein, M. (2015). Global Spatiotemporal Distribution of Soil Respiration Modeled Using a Global Database. Biogeosciences, 12, 4121-4132. https://doi.org/10.5194/bg-12-4121-2015
[8]	Iiyama, I. (2023). Seasonal Change in CO2 Production Rate along Depth in a Grassland Field. Journal of Geoscience and Environment Protection, 11, 106-124. https://doi.org/10.4236/gep.2023.116008
[9]	Iiyama, I., Osawa, K., & Nagata, O. (2012). Soil O2 Profile Affected by Gas Diffusivity and Water Retention in a Drained Peat Layer. Soils and Foundations, 52, 49-58. https://doi.org/10.1016/j.sandf.2012.01.005
[10]	Ise, T., & Moorcroft, P. R. (2006). The Global-Scale Temperature and Moisture Dependencies of Soil Organic Carbon Decomposition: An Analysis Using a Mechanistic Decomposition Model. Biogeochemistry, 80, 217-231. https://doi.org/10.1007/s10533-006-9019-5
[11]	Jassal, R. S., Black, T. A., Novak, M. D., Gaumont-Guay, D., & Nesic, Z. (2008). Effect of Soil Water Stress on Soil Respiration and Its Temperature Sensitivity in an 18-Year-Old Temperate Douglas-Fir Stand. Global Change Biology, 14, 1-14. https://doi.org/10.1111/j.1365-2486.2008.01573.x
[12]	Kaminski, T., Knorr, W., Rayner, P. J., & Heimann, M. (2002). Assimilating Atmospheric Data into a Terrestrial Biosphere Model: A Case Study of the Seasonal Cycle. Global Biogeochemical Cycles, 16, Article No. 1066. https://doi.org/10.1029/2001GB001463
[13]	Kim, G. S., Joo, S. J., & Lee, C. S. (2020). Seasonal Variation of Soil Respiration in the Mongolian Oak (Quercus mongolica Fisch, Ex Ledeb.) Forests at the Cool Temperate Zone in Korea. Forests, 11, Article No. 984. https://doi.org/10.3390/f11090984
[14]	Knorr, W., & Heimann, M. (1995). Impact of Drought Stress and Other Factors on Seasonal Land Biosphere CO2 Exchange Studies through an Atmospheric Tracer Transport Model. Tellus B: Chemical and Physical Meteorology, 47, 471-489. https://doi.org/10.3402/tellusb.v47i4.16062
[15]	Kochy, M., Hiederer, R., & Freibauer, A. (2015). Global Distribution of Soil Organic Carbon—Part 1: Masses and Frequency Distributions of SOC Stocks for the Tropics, Permafrost Regions, Wetlands, and the World. Soil, 1, 351-365. https://doi.org/10.5194/soil-1-351-2015
[16]	Lee, D. K., Doolittle, J. J., & Owens, V. N. (2007). Soil Carbon Dioxide Fluxes in Established Switchgrass Land Managed for Biomass Production. Soil Biology and Biochemistry, 39, 178-186. https://doi.org/10.1016/j.soilbio.2006.07.004
[17]	Liu, X., Wan, S., Su, B., Hui, D., & Luo, Y. (2002). Response of Soil CO2 Efflux to Water Manipulation in a Tallgrass Prairie Ecosystem. Plant and Soil, 240, 213-223. https://doi.org/10.1023/A:1015744126533
[18]	Liu, Z., Deng, Z., Davis, S., & Ciais, P. (2023). Monitoring Global Carbon Emissions in 2022. Nature Reviews Earth and Environment, 4, 205-206. https://doi.org/10.1038/s43017-023-00406-z
[19]	Mahecha, M. D., Reichstein, M., Carvalhais, N., Lasslop, G., Lange, H., Seneviratne, S. I., Vargas, R., Ammann, C., Arain, M. A., Cescatti, A., Janssen, I. A., Migliavacca, M., Montagnani, L., & Richardson, A. D. (2010). Global Convergence in the Temperature Sensitivity of Respiration at Ecosystem Level. Science, 329, 838-840. https://doi.org/10.1126/science.1189587
[20]	Makita, T. (1988). Diffusion Coefficient in Gaseous Phase. In Kagakukougaku-Binran (Ed.), The Society of Chemical Engineers, Japan (pp. 100-104). Maruzen, Tokyo. (In Japanese)
[21]	National Astronomical Observatory of Japan (2020). Rika Nenpyo 2021 (Chronological Scientific Tables 2021) (pp. 532-535). Maruzen Publishing Co., Ltd. (In Japanese)
[22]	Osozawa, S., & Hasegawa, S. (1995). Diel and Seasonal Changes in Carbon Dioxide Concentration and Flux in an Andisol. Soil Science, 160, 117-124. https://doi.org/10.1097/00010694-199516020-00005
[23]	Post, W. M., Emanuel, W. R., Zinke, P. J., & Stangenberger, A. G. (1982). Soil Carbon Pools and World Life Zones. Nature, 298, 156-159. https://doi.org/10.1038/298156a0
[24]	Raich, J. W., & Schlesinger, W. H. (1992). The Global Carbon Dioxide Flux in Soil Respiration and Its Relationship to Vegetation and Climate. Tellus B: Chemical and Physical Meteorology, 44, 81-99. https://doi.org/10.3402/tellusb.v44i2.15428
[25]	Risk, D., Kellman, L., & Beltrami, H. (2002). Soil CO2 Production and Surface Flux at Four Climate Observatories in Eastern Canada. Global Biogeochemical Cycles, 16, Article No. 1122. https://doi.org/10.1029/2001GB001831
[26]	Wang, M., Li, X., Wang, S., Wang, G., & Zhang, J. (2018). Patterns and Controls of Temperature Sensitivity of Soil Respiration in a Meadow Steppe of the Songnen Plain, Northeast China. PLOS ONE, 13, e0204053. https://doi.org/10.1371/journal.pone.0204053
[27]	Wang, W., Chen, W., & Wang, S. (2010). Forest Soil Respiration and Its Heterotrophic and Autotrophic Components: Global Patterns and Responses to Temperature and Precipitation. Soil Biology and Biochemistry, 42, 1236-1244. https://doi.org/10.1016/j.soilbio.2010.04.013
[28]	Xu, L., Baldocchi, D. D., & Tang, J. (2004). How Soil Moisture, Rain Pulses, and Growth Alter the Response of Ecosystem Respiration to Temperature. Global Biogeochemical Cycles, 18, GB4002. https://doi.org/10.1029/2004GB002281
[29]	Zhan, Z., Zhang, R., Cescatti, A., Wohlfahrt, G., Buchmann, N., Zhu, J., Chen, G., Moyano, F., Pumpanen, J., Hirano, T., Takagi, K., & Merbold, L. (2017). Effect of Climate Warming on the Annual Terrestrial Net Ecosystem CO2 Exchange Globally in the Boreal and Temperate Regions. Scientific Reports, 7, Article No. 3108. https://doi.org/10.1038/s41598-017-03386-5
[30]	Zhang, L. H., Chen, Y. N., Zhao, R. F., & Li, W. H. (2010). Significance of Temperature and Soil Water Content on Soil Respiration in Three Desert Ecosystems in Northwest China. Journal of Arid Environments, 74, 1200-1211. https://doi.org/10.1016/j.jaridenv.2010.05.031
[31]	Zhang, Z. S., Dong, X. J., Xu, B. X., Chen, Y. L., Zhao, Y., Gao, Y. H., Hu, Y. G., & Huang, L. (2015). Soil Respiration Sensitivities to Water and Temperature in a Revegetated Desert. Journal of Geophysical Research: Biogeosciences, 120, 773-787. https://doi.org/10.1002/2014JG002805
[32]	Zhou, T., Shi, P., Hui, D., & Luo, Y. (2009). Global Pattern of Temperature Sensitivity of Soil Heterotrophic Respiration (Q10) and Its Implications for Carbon-Climate Feedback. Journal of Geophysical Research, 114, G02016. https://doi.org/10.1029/2008JG000850