SOC is a continuum of progressively decomposing organic components ranging from labile compounds that decompose quickly to recalcitrant compounds that decomposes extremely slowly (Lehmann and Kleber, 2015). Following the SCALE model (Yan et al., 2019), here we categorize SOC into three pools with feedbacks between each other—fast pool (C_l), slow pool (C_h), and microbial biomass (C_b) (Porporato et al., 2003). Combining the biogeochemical transformation, soil erosion/deposition (and resultant landscape evolution), and bioturbation by soil fauna, the SOC mass conservation in a soil column is summarized as Yan et al. (2019):

where subscripts 1 and n denote the surface soil layer and the nth layer below-surface, respectively; C=[Cl,Ch,Cb]T is the SOC content, where C_l, C_h, C_b represent the fast (or litter), slow (or humus), and microbial biomass pools, respectively [ML⁻³]; t represents time [T], z represents soil depth [L], ∇ · q_C is the sum of surface SOC fluxes due to erosion from overland flow and diffusion-like surface transport (Yan et al., 2019); ∇ · [D(z)∇C] is the vertical diffusion process resulting from bioturbation by soil fauna, where D(z) is the bioturbation diffusivity, parameterized as D(z)=Dtope−0.1z (Quijano et al., 2013). The biogeochemical transformation term g represents the carbon influx from above- and below-ground plant residues, and carbon efflux due to decomposition throughout the three pools. Their relationships are shown as Porporato et al. (2003):

where I_litter is the litter input from litterfall and root-litter [ML⁻²T⁻¹]; k_rd is the death rate of microbes [T⁻¹]; r_h is referred to as “isohumic” coefficient (Wild, 1988), which is the fraction of decomposing litter that undergoes humification and ranges from 0.15 to 0.35 [−] (Brady and Weil, 1996; O'dorico et al., 2003; Porporato et al., 2003); r_r defines the fraction of decomposed organic C to CO₂ [−] (0 ≤ r_r ≤ 1 − r_h); K_l and K_h are decomposition rate in fast and slow pool, respectively [T⁻¹]. These two rates are regulated by microbial biomass, soil moisture, C/N ratio, and soil temperature. The influence of soil temperature is not addressed in this study but can be added within this framework as needed.

The control of microbial biomass, soil moisture and C/N ratio on the decomposition rates in the fast pool (K_l) and slow pool (K_h) are given below (Porporato et al., 2003):

where k_l and k_h represent the rates of decomposition as simplified terms that encompass different organic components in fast and slow pools, respectively [L³T⁻¹M⁻¹]; φ is a ratio that is from the reduction of the decomposition rate if the immobilization (controlled by nitrogen content) fails to meet the nitrogen demand by the microbes [−]. φ ≈ 1 in agricultural fields where nitrogen supply is usually sufficient from fertilizers; f_d(θ) [−] represents the soil moisture effects on decomposition. The relationship between soil moisture (θ) and the index f_d(θ) is given as Porporato et al. (2003):

The soil moisture is estimated following the Richard's equation for one-dimensional vertical soil moisture dynamics (Yan et al., 2019). The parameters related to the simulation are based on field measurements (Abaci and Papanicolaou, 2009; Wilson et al., 2018). The optimal soil moisture condition is the field capacity (θ_fc) which provide the maximum decomposition rate (f_d(θ) = 1). Drier or wetter conditions will result in a smaller f_d(θ), and hence reduce the decomposition rate.

Summing the three terms of g in Equation (2) gives the total rate of change of SOC due to biogeochemical transformation:

The total HCL through the soil column and the effective decay rate (D_r) for SOC are expressed as:

Here we focus on the two key processes resulting from agricultural activities —accelerated soil transport rate and tillage-induced mechanical-mixing on top soils. In our simulations, the soil transport (the net flux the of incoming subtracting the outgoing soil material in a unit volume. In this study, it represents soil depth change rate driven by erosion or deposition [L/T]) is a combination of diffusive-like process of topographic elevation, driven by topographic gradient, wind erosion, biogenic disturbance, soil creep, and rain drop splash; and advective erosion of the topography driven by land surface overland flow process. Soil transport not only results in soil erosion and deposition but also causes the spatial redistribution of SOC contents on the land surface. The redistribution rate of SOC is assumed to be positively proportional to the rate of soil transport, which is captured by ∇·q_C term in Equation (1) (Yan et al., 2019).

The tillage happens before planting in Spring season, assuming on day 104 each year, and several days variation around this date does not have any significance bearing on the results presented. The mechanical-mixing from tillage creates a uniform SOC content throughout the tillage depth, assumed here as 20 cm. We investigate the impacts of mechanical-mixing separately from the accelerated soil transport, where the latter arises partially from tillage but also from seasonal landcover change, and other anthropogenic activities.

This study focuses on the CO₂ production by microbial decomposition through the soil column. To further translate this to the amount of CO₂ released from the soil to the atmosphere, one has to incorporate the control of the CO₂ gradient between the soil and atmosphere, mineral absorption, dissolution in soil-water, etc., which is beyond the scope of this study. Also, the impact of soil temperature is an important factor in the SOC decomposition process but not included in this study because the study area has very low local relief (with a gradient about 0.0002), and therefore, the soil temperature is expected to exhibit little spatial difference in HCL over the time scale of 150 years. However, the soil temperature influences could be added within this framework in future studies.

We performed our study in a first-order sub-watershed in the Clear Creek Watershed (CCW) in eastern Iowa (Figure 1). The CCW represents most Midwestern watersheds regarding the Mollisols soils, humid-continental climate (hot summers, cold winters, and wet springs with freeze-thaw periods in late fall and early spring), and predominantly agricultural land use. The surface layer soils are loess-derived, which is highly productive but also highly erodible. The mean annual precipitation is about 889 mm with prominent thunderstorms in the summer and snowfall in the winter. The mean annual temperature is 9°C with mean warmest month at 29°C and mean coldest month at −13°C. This area was prairie wetlands and forests until the European Settlement in 1880s since when farming activities eventually converted about 80% of the lands to crops and pastures. After transforming into farmable lands, soil transport (erosion and deposition) was accelerated since then (Kumar et al., 2018).

The forcing includes daily precipitation from stochastic simulation, daily organic carbon input from plant litters, corn-soybean rotation in alternate years, and soil properties (i.e., saturated hydraulic conductivity, soil porosity, maximum infiltration rate, and Van Genuchten parameters). The forcing that drives the physical and biogeochemical processes is spatially uniform due to a relatively small size study domain (0.15 km²). The initial values of this simulation include light detection and ranging (lidar) Digital Elevation Model (DEM), which provides a 2-m resolution as for the land surface grid size; SOC and soil moisture profiles at each grid cell; and soil layer thickness at each grid cell. The vertical profiles of SOC and soil moisture along the soil depth have seven grids (see values in Table 1) and are updated after each computational time step. The control volumes of the vertical grids are dynamic because of the changing elevations due to surface soil transport (Yan et al., 2019). Except for the elevation, all other variables are spatially uniform at an initial stage of the simulation. The key variables and corresponding initial values are provided here (Table 1). The full list of input parameters are summarized in Yan et al. (2019) and are not repeated here for brevity.

Our simulation starts with an initially exponential SOC profile as a function of depth. This profile is computed based on a steady-state between C input from plant residues and decomposition under climate forcing (Yan et al., 2019) because local prairie wetlands were formed and provided a relatively stable SOC condition after thousands of years after the end of the last glacial episode (Yan et al., 2017). This site has previously served as a validation site of the SCALE model (Yan et al., 2019) in characterizing the spatial heterogeneity of the evolution of SOC profiles. The simulated spatial mean soil transport rate after 150 years of anthropogenic influences is ~1 mm/yr (Yan et al., 2019), which has been validated with field measurements (Wilson et al., 2009; Papanicolaou et al., 2015). We implemented the mechanical-mixing at year 100 when the HCL has reached steady-state from soil transport, litter input, meteorologic forcing, and natural bioturbation driven diffusive mixing. The mechanical-mixing, a result of tillage in agriculture practices, mixes the top 20 cm of the soil creating a uniform SOC profiles but conserves the total carbon mass. Our assumption here is that the soil depth change is all from the surface soil transport and not from the soil bottom weathering. The main reason is that the bottom of soils is layered with loess and glacial drifts, and the weathering processes are much slower than the surface transport rate, which is very different compared to mountainous areas (Anders et al., 2018; Yan et al., 2021).

We simulated a 150-year period to understand the long-term dynamics of HCL (Figure 2) due to the influences of accelerated soil transport (and the resultant landscape evolution) and mechanical-mixing. The HCL reached a dynamic steady-state at around 75 years before the mechanical mixing started, and a new dynamic steady-state re-established approximately 40 years after the implementation of the mechanical mixing at year 100. The accelerated transport of solid phase of soil material leads to an increase of HCL by approximately 40% and 57% before and after the mechanical mixing, respectively (from the initial stage of 28.5 to 40.0 g C/m²/yr at year 75 and 44.8 g C/m²/yr at year 140). To characterize the differences arising from soil erosion and deposition across the landscape, four zones were selected (Figure 2a)—Zones A and D represent areas with the most erosion and deposition, respectively; and Zones B and C are areas with moderate erosion and deposition, respectively. The spatial locations corresponding to the four zones are shown in Figure 2b. The temporally variation of the spatial mean values of each of the four zones as well as the whole watershed are shown in Figure 2c. The annual mean values from year 75 till year 100 for each of the four zones are 21.6, 30.7, 54.5, and 59.3 g C/m²/yr; and from year 140 till year 150 are 16.1, 25.5, 68.5, and 72.8 g C/m²/yr for Zone A, B, C, and D, respectively.

Soil respiration rates from field measurement in temperate ecosystems ranges from 10 to 1,500 g C/m²/yr (Bouma and Bryla, 2000; Raich and Tufekcioglu, 2000; Bond-Lamberty and Thomson, 2010a,b). One study measured the annual average soil respiration rate in a corn and soybean field in central Iowa (Tufekcioglu et al., 2001), which is 600–750 g C/m²/yr, and the mean daily soil respiration in corn and soybean ranges about 0.1 to 5 gC/m²/day. Another estimation in corn in WI, USA (von Haden et al., 2019) shows that the annual accumulated HCL is about 550 gC/m²/year. Assuming that the HCL takes approximately half of the total soil respiration (Raich and Potter, 1995; Bond-Lamberty et al., 2004; Bond-Lamberty and Thomson, 2010a; Hashimoto et al., 2015), the HCL is around 300–550 g C/m²/yr in corn-soybean field in the Midwest U.S., and within the total range of 5–750 g C/m²/yr globally. Even though our simulation results are smaller than the two studies in corn fields, we do find sites where the HCL is over 250 g C/m²/yr (see next subsection). The possible reasons that our simulations show smaller values than the measurement would be that (1) the measurements are at a certain location, and possibly in a lower hillslope location, which may tends to be a higher soil respiration spot in a watershed; (2) our litter input from crop residues is around 32 g C/m²/yr, which may be smaller than other study sites; (3) the autotrophic respiration may take a larger portion than 50% to the total soil respiration, and (4) the annual respiration in other studies are predicted based on a growing season then projected into a year, which may overestimation.

Comparing the spatial variation of the HCL in a watershed, the HCL is higher at depositional sites than erosional sites (Figure 2c). Specifically, soil erosion leads to a reduction of HCL, and deposition leads to an increase in HCL. At erosional sites, the removal of surface soils exposes subsoils, which have less available SOC to produce HCL. Further, the newly exposed subsoils tend to have fast weathering rate, which turns the exposed subsoil layers to be unsaturated with sorbed organic carbon compounds. The high capability to store more SOC of the newly exposed subsoils reduces the release of carbon as HCL. At depositional sites, the near-surface SOC that is removed from erosional sites gradually accumulates at depositional zones. The high SOC availability at depositional sites can enhance the HCL; meanwhile, the burial of soils can suppress the decomposition rates (K_l and K_h) in that the microbial biomass and decomposition coefficients (k_l and k_h, Equation 3) both decrease with depth. The suppression of decomposition would offset the positive influence of the SOC accumulation on the HCL. The simulation result suggests that the accumulation of near-surface SOC in depositional sites plays a more important role on the HCL than the effect of burying SOC. Overall, erosion causes a reduction of HCL by the removal of SOC, which in turn can compensate for the loss of eroded SOC; and deposition leads to an increase of HCL, which offsets the accumulation of deposited SOC.

At year 100 from the start of the simulation (Figure 2c), a mechanical-mixing was implemented to capture the effect of tillage, which vertically homogenized the SOC content of the topsoils (i.e., 20 cm, once a year). The mechanical-mixing reduces the HCL at erosional zones, which eventually slows down the SOC stock loss from erosion. Similarly, the mechanical-mixing enhances the HCL at depositional zones, which eventually slows down the SOC stock accumulation from deposition. In other words, the mechanical-mixing reduces the spatial differences of the rate of SOC stock change among erosional and depositional areas. More discussions about the effect of mechanical-mixing on HCL are given in section 3.4.

Comparing the relationship between the soil transport rate and annual accumulated HCL (Δ Soil Depth) at each 2-D grid box in the study domain, we find that the sites with fast soil deposition rate does not necessarily correspond to high HCL (Figure 3A1), even though depositional sites show higher HCL than erosional sites by average. For example, in Zone D (Figure 3A1, highlighted in blue color) where soil thickness increases most due to deposition, it accumulates more solid soils than other zones, but the burial of SOC suppresses its decay rate (due to smaller decomposition coefficients (k_l and k_h, Equation 3) with increasing soil depth in the model). The highest HCL corresponds to relatively moderate erosional or depositional areas (Zones B and C in Figure 3A1) where there exists steep slopes that are associated with high soil influx (the incoming soil material to a unit volume due to deposition [L/T]) and efflux exchange (the outgoing soil material from a unit volume due to erosion [L/T]).

Specifically, among Zones B and C, even though the soil transport rate (Δ Soil Depth, defined as the net flux of influx minus efflux), is relatively small; locations with high influx and efflux of soil can create an environment that favors high HCL.

In some other areas among Zone B and C where both the influx and efflux of soil are small, the corresponding HCL is also small. We hypothesize that the high variation of influx and efflux among Zone B and C is the reason that causes high HCL variation.

The total HCL has a strong positive monotonic correlation with the SOC redistribution (R² = 0.88, Figure 3B1). At sites where SOC redistribution is positive, meaning a net deposition of SOC, the total HCL increases almost linearly with SOC redistribution. The ratio of total HCL to SOC redistribution is close to but larger than one (Figure 3B1), and the SOC redistribution accounts for approximately 65% of the total HCL in depositional sites. This implies that most of the HCL is from the SOC redistribution, and a small portion of HCL is from local immobile SOC. As with the HCL, the areas that have highest spatial variation of SOC redistribution rates correspond to the very small soil transport rate (Δ Soil Depth) (Figure 4A). This demonstrates that a zero net flux of soil transport would not necessarily result in a zero net flux of SOC redistribution because the SOC concentration is not spatially uniform. There exist areas, even though relatively small, where the values of influx and efflux of soil are similar (where Δsoil depth ≈ 0) but SOC redistribution and HCL variation are relatively high.

Regarding the non-linear relationship between SOC stock and HCL, our results show that the highest HCL is around the moderate SOC stock zone rather than the high SOC stock zone (Figure 3C1). This is because the SOC stock has a positive linear relationship with the soil transport rate (Figure 4B), so the relationship of HCL v.s. soil transport is similar as the one of HCL v.s. SOC stock. In our initial settings, the SOC stock is spatially uniform. The factors that drive the spatial heterogeneity of SOC stock are the SOC redistribution and biogeochemical turnover. The biogeochemical turnover is the outcomes of spatially uniform carbon input from plant residues and the microbial decomposition of organic carbon (represented as HCL in this study). The monotonic positive relationship between soil transport and SOC stock implies that the redistributed SOC happens to largely contribute to the HCL, particularly for high HCL sites sot that the HCL systematically influences the impact of SOC redistribution on SOC stock. In other words, this indicates that the detached and thereafter redistributed SOC can be efficiently and rapidly decomposed by microbial biomass. This non-linear relationship between SOC stock and HCL is different from a study in a farm land in Iowa, which claims that soil respiration has a positive correlation with SOC stock (Tufekcioglu et al., 2001). However, we find the positive correlation holds for below-surface HCL and SOC stock.

The below-surface HCL (the sharp shift of HCL rates, starting about 5 cm below the surface) has one order of magnitude smaller than the total HCL (see section 3.3 for more information). The below-surface HCL shows a high positive relationship with SOC stock (R² = 0.90, Figure 3C2) and a relatively high positive relationship with soil transport (R² = 0.72, Figure 3A2) because the fast soil deposition leads to high SOC stock (Figure 4B). A high soil transport rate not necessarily leads to a high SOC redistribution, and close to zero soil transport rate corresponds to high variation of SOC redistribution (Figure 4A); therefore, the high values of below-surface HCL corresponds to areas where SOC redistribution are relatively small (Figure 3B2). The below-surface HCL shows opposite correlation with soil transport, SOC redistribution, and SOC stock as compared with total HCL. This indicates that the SOC redistribution mostly affects the surface HCL rather than the processes of the below-surface HCL. For the below-surface HCL, it is mostly controlled by SOC stock. We explore other controling factors on the HCL in the following section.

To explore the spatial heterogeneity of the vertical profiles of HCL, the daily mean values of HCL profiles of the four representative zones (Zone A, B, C, and D) are shown in Figure 5A. For all of the zones, the HCL rates decrease rapidly with depth, which is consistent with other studies (Martel and Paul, 1974; Risk et al., 2002; Rumpel et al., 2002; Hirano et al., 2003; Schöning and Kögel-Knabner, 2006; Maier et al., 2011), in that the near-surface HCL is 1 to 2 orders of magnitude higher than those at deeper depths, which is primarily determined by the microbial biomass. The HCL is the product of SOC content and effective decay rates (Equation 6). Similar to the HCL profiles, the near-surface effective decay rate (D_r, in Equation 7) is 1 to 2 orders of magnitude higher than the one below (Figure 5B). D_r is primarily controlled by soil moisture and microbial biomass (Equation 3). The index f_d(θ) (Equation 4) represents the influence of soil moisture on the SOC decomposition in that very dry soil condition can reduce D_r because of low microbial activity, and very wet soil condition can also reduce D_r because of limited oxygen access. The optimal soil moisture condition that favors D_r is the field capacity. We show the vertical profile of f_d(θ), which increases with depth because deeper layer is closer to the field capacity level (Figure 6A). In other words, the high HCL value at near-surface is primarily determined by D_r rather than SOC content, and D_r is primarily determined by the microbial biomass rather than soil moisture. The surface microbial biomass decreases with increasing soil depth (Figure 6B). Studies show that microbial biomass up to 40 cm can be one to two orders of magnitude higher than layers below (Federle et al., 1986; Sun et al., 2021), which is consistent with our simulation results. The reasons may include higher oxygen in the near surface layers that favors the decomposition and physical protection of SOC from microbial that opposes decomposition.

In comparison to the near-surface layer, below-surface HCL profiles are controlled by different factors. The profiles follow a similar trend to the SOC profiles (Figure 1C). At depositional sites, higher content or “noses” show up at deeper depths; whereas at erosional sites, the profiles are exponentially decreasing with depths (Figure 5A). This suggests that the SOC content, rather than the soil moisture or microbial biomass, has a dominant control on the HCL profiles for deep layers. This finding is consistent with the results discussed in section 3.2; and consistent with a study that claims the SOC stock, rather than soil microbial community composition, should be preferentially considered for the spatial soil respiration variance for subsurface layers (Wei et al., 2015). The below-surface effective decay rate, D_r, (Figure 5B) increases slightly toward deeper layers, which is consistent with the f_d(θ) profiles (Figure 6A) and suggests that soil moisture has stronger influence on D_r rather than microbial biomass. The microbial profiles (Figure 6B) also shows a “nose” at deep layers at depositional sites, which indicates the influence from SOC content. Therefore, the SOC content rather than soil moisture and microbial biomass plays a more important role on the below-surface HCL profiles.

To investigate the impacts of mechanical-mixing, resulting from tillage, on the vertical profiles and spatial heterogeneity of HCL, we compare the values before and after the implementation of mechanical-mixing in the model at year 100 (Figure 2c). Since the mechanical-mixing homogenizes the near-surface SOC contents, here assumed for the top 20 cm, it reforms the SOC profiles by reducing the surface content at erosional sites (Figure 7A, Zone A and B) and slightly increasing it at depositional sites (Figure 7A, Zone C and D). With all other conditions the same (including landcover and erosion rate), vertically mixing the near-surface SOC contents results in the reduction of surface values of HCL at erosional sites but addition at depositional sites (Figure 7B). In other words, mechanical-mixing suppresses the CO₂ release at erosional sites but enhances the CO₂ release at depositional sites.

Because the mechanical-mixing enhances the spatial difference of HCL among erosional and depositional sites, the HCL shows a unimodal probability density function (PDF) before mechanical-mixing (Figure 8A) but a bimodal PDF after the mechanical-mixing (Figure 8B). Moreover, among the two PDF plots, the spatial mean values of the HCL is slightly higher after the implementation of mechanical-mixing from 40.0 g C/m²/yr to 44.8 g C/m²/yr due to increase in the vertical distribution of fresh organic carbon input along the soil profile, which provides a source of energy for soil microbes at deeper layers that could stimulate the loss of older buried SOC (Fontaine et al., 2007).