Seasonal Change in CO2 Production Rate along Depth in a Grassland Field

Ippei Iiyama

School of Agriculture, Utsunomiya University, Mine 350, Utsunomiya, Japan.

DOI: 10.4236/gep.2023.116008                

Abstract

Soil is a large terrestrial carbon pool so that the evaluation and prediction of soil respiration is important for understanding and managing carbon cycling between the pedosphere and the atmosphere. For better understanding about characteristics and mechanisms of soil respiration, this study monitored seasonal behaviors of soil gaseous CO₂ concentration profile with relevant soil physical conditions in a meadow field, and numerically analyzed the monitored data sets to inversely determine time-series of depth distributions of CO₂ production rate in the field by assuming optimum ranges of depth and moisture condition for aerobic respiration of soil fauna and flora. The results of the inverse analyses showed that the depth range of intense CO₂ production resided in top soil layers during summer and moved down into subsoil layers in winter, implying that the depth range of main CO₂ sources can change dynamically with seasons. The surface CO₂ emission rates derived from the inverse analyses fell in the range typically found in the same kind of land use. The evaluated mean residence time of gaseous CO₂ in the study field was around half a day. These findings suggested that the modelling assumptions about soil respiration in this study are effective to probe spatial and temporal behavior of respiratory activity in a soil layer, and it is still important to integrate facts about in-situ CO₂ concentration profiles with soil physical parameters for quantitatively predicting possible behaviors of soil respiration in response to hypothetical changes in atmospheric and soil climates.

Keywords

Gas Diffusion Coefficient, Mass Balance Equation, Soil Temperature, Volcanic Ash Soil, Water Content

1. Introduction

Soil respiration is an ecosystem process through which carbon dioxide (CO₂) is released, and oxygen (O₂) is taken in, by plant roots, soil fauna, or soil microbes. Studies on soil respiration have been prompted by growing attention to global carbon budget with the estimations that soils emit 98 ± 12 [PgCO₂ y⁻¹] in the world (Bond-Lamberty & Thomson, 2010) , which can be several times more amount of carbon than fossil-fuel burning and cement manufacture of 6.3 [PgC yr⁻¹] (Schimel et al., 2001) . Therefore, CO₂ emissions from soil surfaces are worth being evaluated intensively under various land uses and soil conditions.

Soil respiration can be both a strong cause and a clear effect of biological activities in soils. For instance, Rao and Ito (1998) reported that root respiration rates of legumes (pigeon pea, chickpea and groundnut) and cereals (sorghum, pearl millet and maize) showed high correlations with nitrogen uptake activities. Liu and Li (2005) gave soil-drying conditions to spring wheat in its early growing season, and suggested that the resultant reduction in root respiration improved crop production through the restriction of utilizing photosynthates. Shimoda et al. (2009) stated that the large soil biomass contributed to a higher respired carbon loss, and the increase in biomass led to a linear increase in ecosystem respiration. Hao and Jiang (2014) also experimentally showed obvious linear relationships between root biomass and soil respiration rate for each month in half a year of a rape growing season. These studies imply the importance of quantifying soil respiration with soil climates in depth for probing, improving, or maintaining geo-environmental conditions.

However, contrary to the knowledge about sizes and features of surface CO₂ emission, studies about below-ground CO₂ production has not been examined so much, though in-situ soil respiration rates occurring along depth can serve as effective knowledge for grasping soil conditions or planning soil managements. In addition, the studies dealing with depth profiles of CO₂ production rates often assumed steady state conditions when analyzing gaseous CO₂ dynamics in soils (Drewitt et al., 2005; Fierer et al., 2005; Kellman et al., 2015) , presumably because it has been thought that temporal behavior of gaseous CO₂ in a soil is so slow that an instantaneous value of respiration can be explained almost thoroughly by a set of above-ground and below-ground conditions at the moment in question. This may result in some weakness in the knowledge about transitional evaluation of respiration activities in in-situ soil profiles.

For excavating the understandings about characteristics and mechanisms of soil respiration, this study aimed at numerically quantifying depth distributions of soil respiration rate with soil climatic conditions. For this aim, depth profile of gaseous CO₂ production rate was modelled as a function of depth and soil moisture condition, with the assumption that there should be a depth at which soil respiration takes place most actively while soil respiration should be restricted under extremely dried- or wetted-conditions. And the parameters of the model were inversely determined by making a transient gas diffusion-reaction equation best reproduce CO₂ concentration profiles measured in a field so that seasonal behaviors of both the depth range and the intensity of soil CO₂ production can be quantitatively analyzed.

2. Materials and Measurements

2.1. Study Site

The study field was a meadow (350 × 80 m; 36˚29'23"N, 139˚59'14"E) located in the Utsunomiya University Farm in Moka city, Tochigi, Japan. The monitoring period in the field was from May 2018 to March 2019. In this period, Moka city had mean temperature of 14.8[˚C] with the total precipitation of 923.5 [mm] including the highest daily precipitation of 60.5 [mm·d⁻¹] on 27 August (Japan Meteorological Agency, 2018) .

The grasses and forbs that had been planted in this field were orchard grass (Dactylis glomerata L. cv. Natsumidori), tall fescue (Festuca arundinacea Schreb. cv. Southern-cross), hybrid ryegrass (Lolium hybridum Hausskn. cv. Tetrelite II), white clover (Trifolium repens L. cv. California Ladino), and alfalfa (Medicago sativa L. cv. Ceres). A row of larch trees (Larix kaempferi (Lamb.) Carr.) had been planted for more than 30 years in the middle of the meadow as wind-breaks, stretching 320 [m] in north-south direction.

Two study sites were set at 2- and 8-meters away from the larch tree row, called the sites A and B, respectively. This study symbolized the sites A and B as the tree-stands area and the harvesting area in the meadow, respectively. The depth range of analysis was from the soil surface to 100 [cm] in depth at each site. Visual observations of the soil profiles in the two sites showed that rooting depth of a larch tree covered almost all the 100-cm soil layer in the site A, while plant roots found in the site B were of grasses and forbs developing mainly within 25 [cm] from the soil surface. These observations induced the prediction that soil respiration activities are more intense and spread more deeply in the site A than in the site B.

The soil group of the both sites was Andosol, which is the second most abundant soil in Japan, covering over 31% in land area (Kanda et al., 2018) , and is found in 44% and 18% of the grasslands and of all the agricultural lands in the country, respectively (Takata et al., 2009) . The soil layers of the sites consisted of a topsoil layer of about 25 [cm] in thick with an underlying transitional layer of about 30 [cm] in thick, followed by a subsoil layer. The soil textures of all the soil layers were classified as clay loam soils based on the soil-texture classification defined by the International Union of Soil Science (IUSS).

2.2. Soil Gas Sampling and Measurements of CO₂ Concentration in Soil Gas Samples

Soil gas samples were taken once a week from gas sampling tubes (Tackett, 1968; Osozawa & Hasegawa, 1995) . Each sampling tube was composed of an inner aluminum-tube through which soil gas was taken out and an outer polyvinyl chloride (PVC) tube that protected the inner tube. The tubes were buried permanently at the two sites for the sampling depths of 10, 20, 40, 50, 90 and 100 [cm]. The size of a soil gas sample was 4 [mL]. A gas sample was enveloped into a disposable plastic syringe (SS-05SZ; Terumo Corporation; Tokyo, Japan) with a three-way stopcock (TS-TR1K; Terumo Corporation; Tokyo, Japan), and was analyzed to determine CO₂ concentration c [Mg·m⁻³] within one hour from the time of sampling.

The measurement system developed by Iiyama & Iimura (2014) was used for determining CO₂ concentration of a soil gas sample. The system was made of a CO₂ probe (GMP343; Vaisala; Vantaa Oyj, Finland), an N₂ gas cylinder for supplying a carrier gas, a digital manometer (PG-100-102VP; NIDEC COPAL ELECTRONICS Corp.; Tokyo, Japan) for monitoring whether the carrier gas flowed steadily, and a personal computer to store the time-series outputs from the CO₂ probe. Since the raw outputs from the CO₂ probe were signals in RS232C protocol, they were transmitted into the personal computer by using a communication software (Excel Logger; C2D2 Co. Ltd.; Nagoya, Japan). A calibration curve for converting the raw data into the values of CO₂ concentration was obtained by using CO₂ standards of 1480, 4980, and 50,200 ppm every time when a set of soil gas samples was to be analyzed.

2.3. Soil Three Phases

Volumetric water content θ [m³·m⁻³] was monitored by using a capacitance-type soil moisture sensing system (Diviner 2000, Sentek Pty Ltd., Stepney, South Australia). The raw data of this system were sampled two to three times a week from the same depths as those for the soil gas samplings, and were converted into the values of θ by using the calibration curve which was obtained specifically for this study field. The details about the specification of this system and the site-specific calibration were described in Iiyama (2016) .

Air-filled porosity a [m³·m⁻³] was evaluated by using the following relation:

$a=1-\theta -{\rho }_d/{\rho }_s$ (1)

where ρ_d [Mg·m⁻³] and ρ_s [Mg·m⁻³] are soil bulk density and soil particle density, respectively. Soil bulk density was measured by weighting the undisturbed soil cores of 50 [cm³] before and after drying at 105[˚C] for 24 [h]. Soil particle density was measured by using pycnometers of Gay-Lussac type and in accord with the Japan industrial standards (JIS) A 1202. These two kinds of density were measured for the depth ranges of 0 - 25, 25 - 55, and 55 - 100 [cm] with triplicate samples, and assumed to have been constant during the study period.

2.4. Soil Gas Diffusion Coefficient

Soil gas diffusion coefficient D_s [cm²·s⁻¹] was calculated as follows:

$D_s=\xi D_a$ (2)

where ξ [non-dim] is soil gas diffusivity and D_a [cm²·s⁻¹] is the diffusion coefficient of gaseous CO₂ in the atmosphere. D_a can be a function of pressure p [kPa] and temperature T_s [˚C] of a location in a soil layer as follows (Campbell & Norman, 1998) :

$D_a\left(p,T\right)=0.139{\left(\frac{T_s+273.15}{273.15}\right)}^{1.75}\frac{p_a}{p}$ (3)

where p_a [kPa] is atmospheric pressure. As for the p [kPa], it was assumed that pressure of soil air is commonly in equilibrium with the atmospheric pressure so that p = p_a anywhere in the soil layer. The measured value of soil temperature was substituted into T_s for each time and place.

The procedure of determining the values of ξ followed the methods of Taylor (1949) , Currie (1960) and Osozawa (1987) , in which an N₂-air binary diffusion process is applied to an undisturbed soil core sample set on a single chamber apparatus. The tracer gas for probing N₂-air binary diffusion process inside the measurement system was O₂, and the varying O₂ concentration in the diffusion chamber was measured by using a Galvanic O₂ sensor (O2-204G; Gastec Corp.; Ayase-City, Japan). The values of ξ were evaluated for a series of air-filled porosities, a, by controlling water content of an undisturbed soil core sample step by step. The sampling depths of the undisturbed soil cores were 10 - 20, 40 - 50, and 90 - 100 [cm] of each site, which were represented the depth ranges of 0 - 25, 25 - 55, and 55 - 100 [cm], respectively, and duplicate samples were taken from each sampling depths.

The measured ξ—a relations were described by a modified expression of the Millington-Quirk model (Millington, 1959; Millington & Quirk, 1961) . The modified expression was as follows:

$\xi \left(a\right)=\{\begin{array}{l}0\left(a<a_0\right)\\ \lambda {\left(a-a_0\right)}^m/{\theta }_s^n\left(a\ge a_0\right)\end{array}$ (4)

where a₀ [m³·m⁻³] is the air-filled porosity ineffective to gaseous diffusion so that ξ takes 0 for a < a₀, and θ_s [m³·m⁻³] is the saturated volumetric water content. The enhancement factor λ [non-dim] and the exponents m [non-dim] and n [non-dim] were fitting parameters for making this model reproduce a measured ξ(a) curve.

2.5. Soil Temperature

Because the evaluations of both K_H and D_s require a value of temperature for each depth and date of soil gas sampling, soil temperatures, T_s [˚C], were monitored during the study period with 30-minute intervals by using micro-loggers (HOBO 8K pendant temperature/alarm (waterproof) data logger (UA-001-08); Onset Computer Corp; Bourne, MA, USA) buried at depths of 5, 10, 20, 40, 70 and 100 [cm]. Then, the time-series of daily mean temperature values were determined by averaging the measured data sets for every date.

3. Numerical Processes

3.1. Governing Equation to Determine Depth Profiles of CO₂ Production Rate

Depth profiles of CO₂ production rates were determined by inversely solving the one-dimensional gaseous diffusion-reaction equation such as:

$\frac{\partial }{\partial t}\left(\left(K_H\theta +a\right)c\right)=-\frac{\partial }{\partial z}\left(-D_s\frac{\partial c}{\partial z}\right)+S$ (5)

where t [s] is time variable and z [cm] is upward-positive vertical location. The surface of the soil layer was set at z = 0, and the spatial domain of analysis was defined as −100 ≤ z ≤ 0. c [Mg·m⁻³] is CO₂ concentration of soil gas as a function of z and t. K_H [non-dim], θ [m³·m⁻³], a [m³·m⁻³], D_s [cm²·s⁻¹], and S [Mg·m⁻³·s⁻¹] are Henry’s constant for gaseous CO₂, volumetric water content, air-filled porosity, soil gas diffusion coefficient, and CO₂ production rate, respectively. The measurements and the formulations about these physical properties will be explained in the following sections.

A problem of solving this partial differential equation is often referred to as a forward problem, since the initial and boundary conditions, the physical properties, and the dependent variable can be regarded as causes, factors, and an effect of a phenomenon of interest, and the procedure to derive c (z, t) goes forward through a from-cause-to-effect process. A depth profile of CO₂ production rate can be determined inversely by making use of the process of solving this forward problem. The inverse analysis in this study was the determination of the term S by optimizing the parameters in the expression of S.

3.2. Henry’s Constant for Gaseous CO₂

Henry’s constant for gaseous CO₂, K_H, takes the values of 1.71, 0.88, 0.53, and 0.36 for the temperature values of 0, 20, 40, and 60 [˚C], respectively (National Astronomical Observatory of Japan, 2020) . Therefore, this study described it as a third-order polynomial function of soil temperature T_s [˚C] so that the polynomial curve best represents the T_s - K_H relations cited above. The resultant expression of K_H was as follows:

$K_H=k_0+k_1{T}_s+k_2{T}_s^2+k_3{T}_s^3$ (6)

where k₀ = 1.71, k₁ = −5.85 × 10⁻², k₂ = 9.75 × 10⁻⁴, and k₃ = −6.25 × 10⁻⁶.

3.3. Model of CO₂ Production Rate

The CO₂ production rate, S [Mg·m⁻³·s⁻¹], represents the activity of aerobic respiration in the soil layer. Since soil moisture condition is one of regulating factors of respiration in a soil, S was formulated as a function of θ. Conceptually, respiration in a soil can be weak with very high soil moisture levels, since aerobic respiration should be depressed with shortage of supply of free oxygen, while very dry conditions should also make soil respiration less active where substrates or enzymes for respiration cannot smoothly travel through soil or be stored sufficiently. Some field studies of soil respiration showed that soil moisture limits respiration only in the lowest and highest extremities of moisture range (Liu et al., 2002; Xu et al., 2004; Jassal et al., 2008; Zhang et al., 2010) . And these past results were supported by other laboratory studies (Guntinas et al., 2013; Zhang et al., 2015) , suggesting that an intermediate and wide range of soil moisture content gives an optimum level of soil respiration whereas the level of it decreases sharply toward both extremities in the soil moisture content range. Because many of the values of soil gas diffusion coefficient in the study sites took practically zero in the first 10% of air-filled porosity, while the most severely dried conditions observed in the study sites occurred with the volumetric water content of around 20%, this study assumed that the respiration activity in the soil layer sharply lowers around the top and bottom 10% to 20% of the range of θ and becomes inactive at both extremities of the θ range, while it remains almost optimum in a broad intermediate range of θ. Based on this assumption, the response of soil respiration to soil moisture condition was modeled with such a function as follows:

$S=\gamma {\rho }_d{s}_{opt}$ (7)

where γ [non-dim] expresses moisture dependency of S while s_opt [Mg·Mg⁻¹·s⁻¹] is the CO₂ production rate per unit mass of dry soil under an optimum soil moisture condition. γ was modeled as follows:

$\gamma =\frac{a}{a+K_g}\frac{\theta }{\theta +K_g}{\left(\frac{{\theta }_s}{{\theta }_s+2K_g}\right)}^{-2}$ (8)

where K_g [m³·m⁻³] is a shape factor of γ. This formulation allows γ to take 1 as its maximum and to take 0 at both θ = θ_s and θ = 0. This study assigned 0.1 to K_g so that γ drops steeply in the first and the last 10% to 20% ranges of θ.

The s_opt was modelled as a depth-dependent value with the assumption that there is a depth at which aerobic microbes and plant roots produce CO₂ most actively. As a first approximation, this feature can be expressed by such a function as follows:

$s_{opt}=\{\begin{array}{l}0\left(z_u\le z\le 0\text{or}z\le z_l\right)\\ s_{max}\frac{z_u-z}{z_u-z_{max}}\left(z_{max}\le z\le z_u\right)\\ s_{max}\frac{z-z_l}{z_{max}-z_l}\left(z_l\le z\le z_{max}\right)\end{array}$ (9)

where z_u [cm] and z_l [cm] are the upper and the lower limits of the depth range of soil respiration, and s_max [Mg·Mg⁻¹·s⁻¹] gives the peak of s_opt at z = z_max. In this study, z_u = 0 [cm] and z_l = −100 [cm] were assigned with the assumption that some respiration activity may be found everywhere in the domain of analysis.

3.4. Inverse Analysis to Determine a Depth Profile of CO₂ Production Rate

A depth profile of CO₂ production rate was evaluated by inversely solving the governing Equation (5). When a problem of solving the governing equation is referred to as a forward problem, and a measured data set of dependent variables, CO₂ concentration c (z, t), of the forward problem is given, such an inverse problem can be defined that some of the parameters in the forward problem would be identified so as to make the outputs of the forward problem best fit the measured data set of c (z, t). As a depth profile of CO₂ production rate was represented by the term S in Equation (5), s_max and z_max were selected as surrogate parameters for the term S, and were adjusted.

Firstly, a subroutine to solve the forward problem (SFP) was made so that it outputs a time-series of depth profiles of c (z, t) for a given set of parameters, initial condition, and boundary conditions. This SFP can be expressed such an input-output system as follows:

$c_{calc}\left(z,t\right)=\text{SFP}\left(cndtns,prmtrs\right)$ (10)

where c_calc (z, t) denotes a time-series of CO₂ concentration profile as a result of calculation in the SFP for given conditions cndtns and parameter values prmtrs. In this context, s_max and z_max were two of the components of the prmtrs.

The main body of the SFP involved a finite difference scheme for solving equation (5). The spatial domain of analysis from 100 [cm] in depth to the soil surface (−100 ≤ z [cm] ≤ 0) was discretized with the increment Δz = 5 [cm]. The temporal domain of analysis was set from t = 0 [s] to a certain upper boundary of the temporal domain, t = t_max [s] with time increment Δt = 1800 [s]. To formulate t = 0 and t = t_max, let c_meas indicate a set of values of c (z, t) obtained from actual measurements. Then, two sets of measured CO₂ concentrations obtained in any two consecutive sampling dates t = t₀ and t = t₁ can be denoted as c_meas (z, t₀) and c_meas (z, t₁). Since soil gas samples were taken once a week, the difference between t₀ and t₁ is equivalent to 7 days. On these definitions about t₀ and t₁, this study assigned the value of t₁ - t₀ to t_max while c_meas (z, t₀) was assigned to the initial condition of the SFP, c_calc (z, 0). A value of c_meas (z', t₀) at any z' without measurement was linearly interpolated by using c_meas (z_a, t₀) and c_meas (z_b, t₀) with such two adjacent vertical locations of measurement z_a and z_b as z_a ≤ z' < z_b. A value of c_meas (z', t₁) was also evaluated in the same manner as c_meas(z', t₀). These arrangements can be expressed as follows:

$c_{calc}\left(z,0\right)=c_{meas}\left(z,t_0\right)$ (11)

$c_{calc}\left(z,t_{max}\right)=c_{calc}\left(z,t_1-t_0\right)$ (12)

The expression (11) gives an initial condition to an inverse analysis. In the expression (12), c_calc (z, t_max) may or may not resemble c_meas (z, t₁), which depends on what combination of values to be assigned to both s_max and z_max in the process (10).

The boundary conditions in the SFP were set as follows:

$c\left(0,t\right)=c_0\left(z=0\left[\text{cm}\right]\right)$ (13)

$\frac{\partial c}{\partial z}=0\left(z=-100\left[\text{cm}\right]\right)$ (14)

where c₀ represents the CO₂ concentration of the atmosphere. This study assigned a value equivalent to 0.04% to c₀ as a practical reference. The lower boundary condition meant a zero-flux boundary condition for imitating the observed behaviors of ∂c/∂z between z = −100 [cm] and z = −90 [cm] in the study sites.

The optimization of s_max and z_max can be formulated as the procedure to search such a combination of s_max and z_max that the difference between c_calc (z, t_max) and c_meas (z, t₁) is minimized. In this study, the difference between c_calc and c_meas was evaluated by applying the concept of the least square:

$Err=\sqrt{\underset{i=1}{\overset{6}{{\displaystyle \sum }}}{\left(c_{calc}\left(z_i,t_{max}\right)-c_{meas}\left(z_i,t_1\right)\right)}^2}$ (15)

where Err is an error criterion about the difference between a calculated and a measured depth profiles of CO₂ concentration. i is an index for a vertical location of measurement, by which the six depths of soil gas sampling are differentiated. Then, the process to optimize s_max and z_max can be recognized as the process to find a certain pair of s_max and z_max that minimizes Err in Equation (15) in the following steps:

Step 1. Set up the data sets of c_meas, cndtns, and prmtrs except for s_max and z_max. As components of cndtns, values of an air-filled porosity and a soil temperature for a time and a location were evaluated by linearly interpolating the measured values obtained through a set two consecutive measurements at two neighboring measurement depths.

Step 2. Define domains of s_max and z_max for searching an optimum pair “ $s_{max}^{\ast }$ and $z_{max}^{\ast }$ ” that minimizes Err.

Step 3. Select a tentative pair of s_max and z_max from the domains defined in step 2.

Step 4. Call SFP to generate c_calc by using the data sets given in steps 1 and 3.

Step 5. Evaluate Err.

Step 6. If the Err evaluated in Step 5 is smaller than the minimum of ever recorded, then update a best-ever pair $s_{max}^{\ast }$ and $z_{max}^{\ast }$ with the currently-tried s_max and z_max, and the minimum-ever Err is also updated with the Err just evaluated in Step 5. If not, do nothing.

Step 7. If all the candidates of s_max and z_max in the domains defined in step 2 have been tried, exit the repetition of the process between step 3 and step 7. If not, go back to step 3.

At the end of the process mentioned above, the optimum pair of s_max and z_max has been obtained as $s_{max}^{\ast }$ and $z_{max}^{\ast }$ . By applying this searching scheme to every time-interval of measurement of CO₂ concentration profile in the field, a time-series of depth profile of S was determined.

After the series of inverse analyses were implemented, time-series of CO₂ emission rates from the soil surfaces in the sites A and B were also evaluated by calculating the series of c_calc (z, t) with optimum pairs of s_max and z_max for all the time-intervals of measurement of CO₂ concentration profile in the field. A surface CO₂ emission rate q_sur for each time-interval of measurement of CO₂ concentration profile was, then, evaluated by the following expression:

$q_{sur}=\frac{1}{J}{\displaystyle {\sum }_{j=0}^J\left(-D_s\frac{c_{calc}\left(z=0\left[\text{cm}\right],t_j\right)-c_{calc}\left(z=-5\left[\text{cm}\right],t_j\right)}{0-\left(-5\right)\left[\text{cm}\right]}\right)}$ (16)

where j denotes an index of time increment for the finite difference scheme for solving Equation (5), J is the maximum of j such that t_J - t₀ equals to a time-interval of measurement of CO₂ concentration profile, being 7 days in this study. c_calc (0, t_j) and c_calc (−5, t_j) are the CO₂ concentrations calculated for a time t_j and the vertical locations of z = 0 [cm] and z = −5 [cm], respectively. D_s was evaluated in the same manner as stated in the section 2.4.

4. Results

4.1. Time-Series of CO₂ Concentration Profiles

The time-series of CO₂ concentrations at each gas sampling depth are shown in Figure 1. Through the entire study period, the values of CO₂ concentration in both sites increased with depth particularly in the top soil layers, while the gradients of CO₂ concentration were less steep in the subsoil layers and there was almost no difference in CO₂ concentration between 90 and 100 [cm] in depth. These observed facts suggested that gaseous CO₂ were produced mainly in the top soil layers in the warm seasons while there was almost no CO₂ gas exchange between the soil layers below and above 100 [cm] in depth.

The CO₂ concentrations in both sites had increased successively by the end of July when the values at 100 [cm] in depth exceeded 3% in the site A and 2.5% in the site B. The CO₂ concentrations decreased in August, and increased again toward their annual maximum values that came out in September. From October to February, the CO₂ concentrations in the two sites had monotonously decreased to 1% or less. In March, only in the site A, there was a sign of recovery in CO₂ concentration at each depth, presumably because the end of dormancy of the larch trees was earlier than the spring revival of the grasses and forbs.

4.2. Soil Three Phases

The time-series of soil water contents are plotted in Figure 1(iii-A) and Figure 1(iii-B) for the sites A and B, respectively. The measurement depths in the top soil layers in both sites had experienced substantial wetting and drying cycles particularly in the first half of the study period. In these measurement depths, the soil water contents gradually decreased day by day, and steeply recovered every time when there was a rainfall of 20 [mm·d⁻¹] or more. In August, the water contents in the top soil layers dropped to the lowest level of the year. Then, the rainfall of 60.5 [mm·d⁻¹] on 8/27 clearly recharged soil water storage in each site, followed by gradual increase in soil water content with continual rainfalls. In contrast to the top soil layers, the water contents in the subsoil layers were relatively stable, and seldom reduced under 0.5. From autumn to the latter days in winter, the depth distributions of soil water content were stable in both sites, because of exceptionally small amount of precipitation with low evaporation demands in those seasons.

The soil bulk densities and particle densities of the two sites are listed on Table 1. By inputting these values and the measured soil water contents into Equation (1), air-filled porosities were determined and used for calculating Equations (4), (5), and (8).

Figure 1. The time-series data sets of (i) daily precipitation, (ii-A) and (ii-B) CO₂ concentration, (iii-A) and (iii-B) volumetric water content, and (iv-A) and (iv-B) soil temperature. Either of “A” and “B” in a sub-index of a sub-graph label denotes that the chart is for the data set from either of the sites A and B.

Table 1. The soil bulk densities and particle densities of the two sites. Each value is an average for triplicate samples with a standard deviation in the parentheses.

4.3. Soil Temperature

Figure 1(iv-A) and Figure 1(iv-B) show the time-series of daily soil temperatures in the sites A and B, respectively. The temperature regimes in the two sites were like each other. The whole soil layers had almost continuously warmed up until the top most soil layers reached more than 27[˚C] in late July. Then, the shallower layers started cooling down while the deeper layers continued warming, and in the late September, the former vertically ascending temperature profile turned into the temperature profile with vertically descending order. After the reversal of the temperature gradient, the subsoil layers were consistently warmer than the top soil layers by the end of January and, then, the temperature profiles gradually varied toward vertically uniform ones in February and March to go to the next reversal of the temperature gradient after the study period.

4.4. Soil Gas Diffusion Coefficient

Figure 2 describes the soil gas diffusivities in relation to air-filled porosity for the undisturbed soil cores sampled in the study field. The maximum of all the measured gas diffusivities was merely 0.08, suggesting that both sites did not contain so many easy pathways for soil gas species as expected from their high total porosities mainly because most of the soils in this study usually had water contents of more than 0.5. Under very wetted conditions, the soil gas diffusivities in both sites got close to zero in the first 10% of air-filled porosity, and, with a few exceptions, were completely zero with air-filled porosity being less than 5%, suggesting that some portions of soil air-phase were hardly accessible to the atmosphere. The measurement results also indicated that the gas diffusivities of the top soil layers were higher than those of the subsoil layers for a given air-filled

Figure 2. The soil gas diffusivities in relation to air-filled porosity for the undisturbed soil cores sampled in the study field. Three sub-graphs (i), (ii), and (iii) in the upper row describe the data sets from the site A, while the lower three sub-graphs (iv), (v), and (vi) are for the site B. Each sub-graph is labelled with a depth range in which a data set represents the gas diffusivities. A set of values of the parameters in equation (4) written on each sub-graph are obtained by the least-square fitting of equation (4) to the measured data set on the sub-graph, and are used for drawing the solid line in the sub-graph.

porosity, as the soils in the site A often got dry comparatively easily to those in the site B, owing to larger amounts of macro-pores developed by tree standings (Iiyama & Hirai, 2014) .

The modelled curves of Equation (4) expressed each of the measured data plots well, implying that the mathematical expression of Equation (4) was suitable for reproducing such features of the soils in this study that high ordinary water content prevented soil gas from easily travelling through the soil matrix even though the soil had very high porosity. The parameter sets of the modelled curves in Figure 2 served as inputs for the inverse analyses of this study.

4.5. Inversely-Determined CO₂ Production Rates

Figure 3(i-A) and Figure 3(i-B) show the time-series of the depth distributions of CO₂ production rate in the sites A and B, respectively. The physical unit of the values of CO₂ production rate in the chart is converted to grams of CO₂ produced per square meter of horizontal plane per 5-cm thickness of soil layer per day, “gCO₂ m⁻² (5 cm)⁻¹·d⁻¹”.

Overall, the calculation results reflected phenological traits commonly found on both arboreal and herbaceous plants in the study field. For example, from May to July, the CO₂ production rates in the top soil layers in both sites had clearly increased with soil temperature rising there during that period. And after the end of September, the intensity of CO₂ production at each site became weaker as the soil temperature in the top soil layers were lowering.

5. Discussion

5.1. Seasonal Behaviors of CO₂ Production Rate Profiles

The depth range with the highest level of CO₂ production rate had been found in the top 20 [cm] of the soil layers in both sites by the end of July, suggesting that plant roots and soil microbes had acted the most lively in the 20 [cm] layer during warming seasons. Then, a sudden reduction in CO₂ production rates occurred coincidently in the two sites in August. This was likely because there was no rain in the study field for more than half a month so that soil water content in the top soil layer fell to the lowest of the year, having been insufficient to maintain an ordinary level of respiratory activity. After the tentative reduction of CO₂ production, the most active CO₂ sources again emerged in the top 20 cm in depth and stayed there until September. From October, the location with the highest value of CO₂ production rate became deeper at each site as the temperature gradients along depth reversed before the beginning of October, and the subsoil layers became warmer than the top soil layers in the latter half of the whole study period.

Some site-specific features also appeared. The results of calculation about the site A showed that the CO₂ production rates near the tree stands had dropped down in the autumn and stayed in the lowest level in the winter. In the study field, autumnal tints began on the larch trees from around 11/20, and almost all their leaves had fallen by 12/15, meaning that the period through which the evaluated soil respiration was weakening coincided with the progress of coloring the larch tree leaves. This coincidence implied that the temporal change in CO₂ production rate in the soil layer can explain the transition in stages of development in plants.

Figure 3. The time-series of (i) the depth distributions of CO₂ production rate, (ii) the surface CO₂ emission rates, and (iii) the amount of gaseous CO₂ stored in the soil layers of 100 [cm] in thick in the sites A and B. The data sets on the sub-graphs (i)-A and (i)-B are for the sites A and B, respectively, and the physical unit of the numerals in the legend of each sub-graph is grams of CO₂ produced per square meter of horizontal plane per 5-cm thickness of soil layer per day, “gCO₂ m⁻² (5 cm)⁻¹·d⁻¹”.

On the other hand, the results of calculation about the site B indicated that the CO₂ production in the harvesting area had been weaker than those near the tree stands, obviously reflecting the shallower and narrower elongation of roots of the grasses and forbs than that of the tree stands. For example, in the site B, the depth range with more than 1 [gCO₂ m⁻² (5 cm)⁻¹·d⁻¹] did not go beyond 40 cm in depth, while that with similar amounts of CO₂ production in the site A expanded into 80 [cm] beneath the surface at maximum. In addition, this numerical result also allowed to infer that crop harvesting was done before the grasses and forbs started becoming less active, since they were harvested on 5/24, 7/25, and 9/24 in that year, while the evaluated CO₂ production rates in the harvesting area started decreasing after the end of September.

5.2. Surface CO₂ Emission Rate and Gaseous CO₂ Storage

Figure 3(ii) shows the time-series of the surface CO₂ emission rates in the two sites, which were derived by using Equation (16). The surface CO₂ emission rate in the site A had been roughly twice as large as that in the site B at any time in the study period, due mainly to larger size of respiratory activities by the tree stands than the grasses and forbs. In the site A, the maximum of the values in the study period was 35 [gCO₂ m⁻²·d⁻¹], almost equivalent to 35 [MgC ha⁻¹·yr⁻¹], found in the end of July, while that in the site B was 15 [gCO₂ m⁻²·d⁻¹]. And, even during the winter, the CO₂ production continued, resulting in the 11-month averages of 17.1 and 6.5 [gCO₂ m⁻²·d⁻¹] in the sites A and B, respectively. These values were comparable with annual surface CO₂ emissions observed on similar land uses in past studies (Risk et al., 2002; Fierer et al., 2005; Lee et al., 2007) and, thus, were likely to be realistic.

Figure 3(iii) shows the amount of gaseous CO₂ stored in the soil layers of 100 [cm] in thick in the two sites. In the site A, the gaseous CO₂ had increased in the warming seasons up to 15 [gCO₂ m⁻²] while the yearly maximum value found in the site B was 7 [gCO₂ m⁻²]. Averaging the time-series data sets, the amounts of gaseous CO₂ contained in the soil layers were 7.6 [gCO₂ m⁻²] and 4.0 [gCO₂ m⁻²]. Compared these values with the surface CO₂ emission rates evaluated above, the mean residence time of gaseous CO₂ in the study field were estimated as 0.44 [d] and 0.62 [d] in the sites A and B, respectively, suggesting that it is important for managing this type of field to maintain smooth exchange between soil air and the atmosphere for proper growth of plants and soil fauna.

In summation, the results of the inverse analyses of this study suggested that even simple assumptions about soil respiration as adopted in Equations (7), (8), and (9) make it possible to clarify environmental behaviors of seasonal respiratory activities in soils when adequate mass-balance and transport models are applied to time-series data sets of depth distribution of CO₂ concentration monitored in a field. It was also implied that once a set of relevant parameters in a governing equation is confirmed, it will be possible to predict quantitatively and mechanistically spatial and temporal behaviors of soil respiration induced by hypothetical changes in atmospheric and soil climates.

6. Conclusion

This study aimed at quantifying and characterizing soil respiration in the form of the depth distribution of CO₂ production rates. For this aim, the depth distributions of soil gaseous CO₂ concentration with relevant soil physical parameters were monitored in a meadow, and were inversely analyzed by using a one-dimensional soil gas transport-production equation. The soil gas diffusivity was modelled by a modified Millington-Quirk model. The CO₂ production rate was modelled with such two surrogate parameters as the depth at which respiration occurs most actively in a soil layer and the intensity of the most active soil respiration, and the determination of these two parameters was made the solution of the inverse analysis. The findings from the measurements and the analyses were as follows:

• The CO₂ production model proposed in this study performed adequately as a first-approximation of spatial and temporal distribution of aerobic respiratory activities in a soil layer. Seasonal behaviors of the depth profiles of CO₂ production rates were understood quantitatively by being aware of the changes in the two surrogate parameters of the model. For instance, the depth range of intense CO₂ production was found in the top soil layers during summer while it moved down into the subsoil layers in winter, implying that main CO₂ sources can change with seasons from respiratory activities in a top soil layer to those in deeper parts of the soil layer.

• Surface CO₂ emission rates were derived from the inversely determined CO₂ production rates. The resultant values were comparable with the surface CO₂ emissions from the same types of land use reported in past studies, suggesting that the inversely evaluated values were likely to be reasonable.

• The mean residence time of gaseous CO₂ in the study field was also estimated by using the results of measurements and analyses of this study. The resultant values were around half a day, and suggested that degradation in soil gas transportability can easily inhibit aerobic respiration in a soil layer like this study field and, therefore, it is important to evaluate the rate of CO₂ production in parallel with the size of CO₂ pool in a soil for properly managing soil gaseous environment.

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]	Bond-Lamberty, B., & Thomson, A. (2010). Temperature-Associated Increases in the Global Soil Respiration Record. Nature, 464, 579-582. https://doi.org/10.1038/nature08930
[2]	Campbell, G. S., & Norman, J. M. (1998). An Introduction to Environmental Biophysics (2nd ed., p. 279). Springer.
[3]	Currie, J. A. (1960). Gaseous Diffusion in Porous Media Part I—A Non-Steady State Method. British Journal of Applied Physics, 11, 314-317. https://doi.org/10.1088/0508-3443/11/8/302
[4]	Drewitt, G. B., Black, T. A., & Jassal, R. S. (2005). Using Measurements of Soil CO2 Efflux and Concentrations to Infer the Depth Distribution of CO2 Production in a Forest Soil. Canadian Journal of Soil Physics, 85, 213-221. https://doi.org/10.4141/S04-041
[5]	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
[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]	Hao, Q., & Jiang, C. (2014). Contribution of Root Respiration to Soil Respiration in a Rape (Brassica campestris L.) Field in Southwest China. Plant, Soil and Environment, 60, 8-14. https://doi.org/10.17221/425/2013-PSE
[8]	Iiyama, I. (2016). Differences between Field-Monitored and Laboratory Measured Soil Moisture Characteristics. Soil Science and Plant Nutrition, 62, 416-422. https://doi.org/10.1080/00380768.2016.1242367
[9]	Iiyama, I., & Hirai, T. (2014). Subsoil Water Available in Suction and Poor in Amount under Continuous Grassland Use. Soil Science and Plant Nutrition, 60, 439-447. https://doi.org/10.1080/00380768.2014.915199
[10]	Iiyama, I., & Iimura, D. (2014). A Simple CO2 Gas Analyzing System for Soil Gas Samples. Journal of Japanese Society of Soil Physics, 128, 33-38. (In Japanese with English Summary)
[11]	Japan Meteorological Agency (2018). Weather Station Data in Moka City in Utsunomiya Local Meteorological Office. http://www.jma.go.jp/jma/en/quickinfo/quickinfo.html
[12]	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, 1305-1318. https://doi.org/10.1111/j.1365-2486.2008.01573.x
[13]	Kanda, T., Takata, Y., Kohyama, K., Ohkura, T., Maejima, Y., Wakabayashi, S., & Obara, H. (2018). New Soil Maps of Japan Based on the Comprehensive Soil Classification System of Japan—First Approximation and Its Application to the World Reference Base for Soil Resources (2006). Japan Agricultural Research Quarterly, 52, 285-292. https://doi.org/10.6090/jarq.52.285
[14]	Kellman, L., Myette, A., & Beltrami, H. (2015). Depth-Dependent Mineral Soil CO2 Production Processes: Sensitivity to Harvesting-Induced Changes in Soil Climate. PLOS ONE, 10, e0134171. https://doi.org/10.1371/journal.pone.0134171
[15]	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
[16]	Liu, H.-S., & Li, F.-M. (2005). Photosynthesis, Root Respiration, and Grain Yield of Spring Wheat in Response to Surface Soil Drying. Plant Growth Regulation, 45, 149-154. https://doi.org/10.1007/s10725-004-7864-6
[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]	Millington, R. J. (1959). Gas Diffusion in Porous Media. Science, 130, 100-102. https://doi.org/10.1126/science.130.3367.100.b
[19]	Millington, R. J., & Quirk, J. P. (1961). Permeability of Porous Solids. Transactions of the Faraday Society, 57, 1200-1207. https://doi.org/10.1039/tf9615701200
[20]	National Astronomical Observatory of Japan (2020). Rika Nenpyo 2021 (Chronological Scientific Tables 2021) (pp. 532-535). Maruzen Publishing Co. Ltd. (In Japanese)
[21]	Osozawa, S. (1987). Measurement of Soil-Gas Diffusion Coefficient for Soil Diagnosis. Soil Physical Conditions and Plant Growth, 66, 53-60. (In Japanese with English summary)
[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]	Rao, T. P., & Ito, O. (1998). Differences in Root System Morphology and Root Respiration in Relation to Nitrogen Uptake among Six Crop Species. Japan Agricultural Research Quarterly, 32, 97-103.
[24]	Risk, D., Kellman, L., & Beltrami, H. (2002). Soil CO2 Production and Surface Flux at Four Climate Observatories in Eastern Canada. Global Biogeochemical Cycles, 16, 69-1-69-12. https://doi.org/10.1029/2001GB001831
[25]	Schimel, D. S., House, J. I., Hibbard, K. A., Bousquet, P., Ciais, P., Peyllin, P., Braswell, B. H., Apps, M. J., Baker, D., Bondeau, A., Canadell, J., Churkina, G., Cramer, W., Denning, A. S., Field, C. B., Friedlingstein, P., Goodale, C., Heimann, M., Houghton, R. A., Melillo, J. M., Moore, B., III, Murdlyarso, D., Noble, I., Pacala, S. W., Prentice, I. C., Raupack, M. R., Rayner, P. J., Scholes, R. J., Steffen, W. L., & Wirth, C. (2001). Recent Patterns and Mechanisms of Carbon Exchange by Terrestrial Ecosystems. Nature, 414, 169-172. https://doi.org/10.1038/35102500
[26]	Shimoda, S., Lee, G., Yokoyama, T., Liu, J., Saito, M., & Oikawa, T. (2009). Response of Ecosystem CO2 Exchange to Biomass Productivity in a High Yield Grassland. Environmental and Experimental Botany, 65, 425-431. https://doi.org/10.1016/j.envexpbot.2008.12.007
[27]	Tackett, J. L. (1968). Theory and Application of Gas Chromatography in Soil Aeration Research. Soil Science Society of America Proceedings, 32, 346-350. https://doi.org/10.2136/sssaj1968.03615995003200030025x
[28]	Takata, Y., Nakai, M., & Obara, H. (2009). Digital Soil Map of Japanese Croplands in 1992. Japanese Journal of Soil Science and Plant Nutrition, 80, 502-505. (In Japanese)
[29]	Taylor, S. A. (1949). Oxygen Diffusion in Porous Media as a Measure of Soil Aeration. Soil Science Society of America Proceedings, 14, 55-61. https://doi.org/10.2136/sssaj1950.036159950014000C0013x
[30]	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, Article No. GB4002. https://doi.org/10.1029/2004GB002281
[31]	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
[32]	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