Soil Organic Carbon Dynamics in Semi-Arid Irrigated Cropping Systems

Abstract

Insufficient characterization of soil organic carbon (SOC) dynamics in semi-arid climates contributes uncertainty to SOC sequestration estimates. This study estimated changes in SOC (0–30 cm depth) due to variations in manure management, tillage regime, winter cover crop, and crop rotation in southern Idaho (USA). Empirical data were used to drive the Denitrification Decomposition (DNDC) model in a “default” and calibrated capacity and forecast SOC levels until 2050. Empirical data indicates: (i) no effect (p = 0.51) of winter triticale on SOC after 3 years; (ii) SOC accumulation (0.6 ± 0.5 Mg ha^–1 year^–1) under a rotation of corn-barley-alfalfa^x3 and no change (p = 0.905) in a rotation of wheat-potato-barley-sugarbeet; (iii) manure applied annually at rate 1X is not significantly different (p = 0.75) from biennial application at rate 2X; and (iv) no significant effect of manure application timing (p = 0.41, fall vs. spring). The DNDC model simulated empirical SOC and biomass C measurements adequately in a default capacity, yet specific issues were encountered. By 2050, model forecasting suggested: (i) triticale cover resulted in SOC accrual (0.05–0.27 Mg ha^–1 year^–1); (ii) when manure is applied, conventional tillage regimes are favored; and (iii) manure applied treatments accrue SOC suggesting a quadratic relationship (all R² > 0.85 and all p < 0.0001), yet saturation behavior was not realized when extending the simulation to 2100. It is possible that under very large C inputs that C sequestration is favored by DNDC which may influence “NetZero” C initiatives.

1. Introduction

Soil organic carbon (SOC) accumulation is of interest in agroecosystems as a gauge of relative soil quality through benefitting physical soil properties and influencing soil biogeochemistry. Interest in SOC as a pool for global C sequestration continues as the Intergovernmental Panel on Climate Change (IPCC) promotes increasing SOC as part of integrated response options to mitigate global C emissions [1]. The magnitude of mitigation that SOC storage can provide is dependent on climate, soil, and agroecological conditions [2]. For example, there is a positive relationship between SOC content and fine particle size silts and noncrystalline clays [3,4,5]. Some climates are predisposed for SOC sequestration while others are less inclined to sequester C [6]. Complicating this, various management practices such as varied tillage regimes, winter cover during annual cropping rotations, and manure, biosolids, or nitrogen applications have been shown to affect SOC accumulation [7,8,9]. Therefore, the importance of discussing SOC accumulation contextually among various climates and management strategies cannot be understated.

Soil organic C sequestration estimates are not abundant in semi-arid environments. These environments have been estimated to cover ~18% of global land surface area and will contribute unnecessary uncertainty to global SOC sequestration estimates if not properly characterized [10]. Soil C sequestration in semi-arid cropland has been studied in some areas of the United States [11] but estimates for specific regions could be improved. The cold semi-arid region of the United States, Köppen-Geiger climate classified “BsK”, extends from eastern Montana south through eastern Colorado and is also present along the Snake River Plain in southern Idaho [12]. In this region of Idaho, agriculture consists of dryland cropping systems of barley (Hordeum vulgare) and wheat (Triticum aestivum) or more diverse irrigated systems comprised of potato (Solanum tuberosum), sugar beet (Beta vulgaris ssp.), dry beans (Phaseolus vulgaris), corn (Zea mays), and alfalfa (Medicago sativa) production. Studies of SOC dynamics are lacking in this region. A few have considered dryland cropping systems [13,14] or native vegetation [15], however, knowledge of SOC dynamics under irrigated cropping systems in this region can be improved.

Agricultural computer models are used to simulate real world systems in order to obtain inferences where observed data is absent. Several biogeochemistry models commonly utilized by researchers are the Environmental Policy Integrated Climate Model (EPIC), Rothamsted carbon model (RothC), denitrification decomposition model (DNDC), daily century model (DayCent), and the root zone water quality model (RZWQM2) [16,17,18,19,20]. Typically, these system models are used in “calibration-validation” cases in which models are calibrated to a subset of real observations and validated against additional observations to gauge the accuracy of the model in representing a particular system. Model calibration and validation for a particular use case is imperative if accurate forecasting of data is anticipated. However, as these models are openly available to the public they are sometimes being used in an “off-the-shelf” capacity without calibration to make inferences about potential for C storage as part of Net-Zero C neutrality initiatives. To that end, the performance of these models used in such a capacity is relatively unknown. One popular model in use, due to its longevity and graphical user interface, is the DNDC model.

The DNDC model is a process-based biogeochemistry model in which users specify climate parameters (temperature and precipitation), soil properties (slope, texture, bulk density, pH, initial SOC), and management practices (crop, tillage, fertilization and manure, irrigation, grazing, etc.) to drive a simulation. The DNDC model operates daily to calculate crop growth and balance C and N pools from the interaction of individual submodules. Observations can be used to supplement model information as available; required inputs include minimal climate data (daily average temperature and precipitation), soil data (pH, and initial SOC), and management practices (crop grown, planting and harvest dates, tillage, etc.) For a complete description of the DNDC model refer to the user manual [21] and a guide to DNDC model version and development [22].

The opportunity exists to improve estimates of SOC accumulation under irrigated semi-arid environments in southern Idaho due to the presence of diverse research sites, including one within the United States Greenhouse Gas Reduction through Agricultural Carbon Enhancement Network initiative (GRACEnet). Three long-term studies located in southern Idaho are relevant to this work. The Long-Term Manure study, hereafter LT Manure, was designed to identify the effect of manure application rate and timing (annual vs. biennial application) on nutrient cycling and greenhouse gas emissions. The CoverCrop study was designed to determine the effect of integrating cover crops into continuous corn rotations utilizing conventional (disk/chisel) and minimum-till (strip-tillage) practices on nutrient cycling. The third study is part of the GRACEnet initiative, referred to as GRACEnet throughout, to assess diverse fertilizer forms and timing on nutrient cycling and greenhouse gas emissions. These studies provide a diverse set of management practices to improve estimates of SOC dynamics in irrigated semi-arid climates. Additionally, these studies provide a challenging array of systems to assess the performance of biogeochemistry models in a default and calibrated capacity. Therefore, objectives of this work were (i) to estimate SOC changes due to variation in manure management, tillage regime, cover crop, and crop rotation; (ii) to assess the performance of the DNDC model with minimal “default” input and after calibration for SOC, crop biomass C, N mineralization, soil water contents, and actual evapotranspiration (AET); and (iii) to use the calibrated model to make inferences on the long-term impact on SOC of diverse management practices taking place on irrigated croplands in the semi-arid region of southern Idaho.

2. Materials and Methods

2.1. Sites

All studies are located within the greater Twin Falls area in ID (42°15′0” N, 114°30′0” W). Mean annual temperature and precipitation are 9.6 °C and 24 cm year^–1, respectively. Soils of the region can be characterized as loess deposits overlying basalt; typical taxonomy is represented by the Portneuf soil series (coarse-silty, mixed, superactive, mesic Durinodic Xeric Haplocalcids). A recent study reported some basic soil properties of the region (n = 75); mean soil pH was 8.0, silt and clay size particle fraction was 538 and 315 g kg^–1 respectively, soil organic carbon content was 13 mg kg⁻¹, and soil carbonate content was 59 g kg^–1 [23].

The LT manure study was initiated in the fall of 2012 with a 4 year crop rotation of wheat–potato–barley–sugarbeet. Plots were 18.3 m by 12.2 m and arranged in a randomized complete block design with four replications. Treatments include a no-treatment control (Control), spring application of synthetic fertilizer (Fertilizer), and drystack dairy manure applied at rates of 18, 36, or 52 Mg ha^–1 on a dry weight basis annually or biennially (18A, 18B, 36A, 36B, 52A, and 52B, respectively) in the fall. Synthetic fertilizer applications were made based on pre-plant soil analysis following University of Idaho guidelines for each crop as well as in-season petiole sampling for potato [24]. All plots were disked following manure and synthetic fertilizer applications to incorporate treatments; plots were moldboard plowed preceding potato and sugarbeet crops. As nutrient requirements were intended to be met each year, synthetic fertilizer applications were made in some manure treatments in some years. Additional information regarding study set-up is reported in Leytem et al. [25].

The CoverCrop study was initiated in the fall of 2015 as a continuous corn cropping system. Plots 12.2 m by 12. 2 m were established using a two by four split plot design with a main treatment of tillage (disk/chisel plow and no-till direct seeding/strip tillage) and a secondary treatment of cover crop by manure application (winter triticale with manure, winter triticale without manure, fallow with manure, fallow without manure) with four replications. For the purposes of this study, disk/chisel plow use was considered conventional tillage and strip tillage was considered minimum tillage. Hereafter, treatments will be referred to by hyphenated character code; conventional tillage (CT) or minimum-till (MT), manure application (M) or no manure application (NM), and triticale (T) or fallow (F). Drystack dairy manure was applied in manured treatments each fall at a target weight of 52 Mg ha^–1 on dry weight basis before cover crop seeding. Treatments without manure application received synthetic fertilizer in the spring based on pre-plant soil testing and University of Idaho guidelines for corn [24]. Where manure applications did not meet nutrient recommendations, synthetic fertilizer was supplemented accordingly. All conventionally tilled plots were disked following manure or synthetic fertilizer application to incorporate treatments.

The GRACEnet study was initiated in the fall of 2012 using a typical regional dairy forage rotation of corn-barley-alfalfa-alfalfa-alfalfa. Plots 21.3 m by 22.9 m were arranged in a randomized complete block design with four replications of six treatments; fall applied dairy manure (Fall manure), fall applied composted dairy manure (Fall compost), spring applied dairy manure (Spr manure), spring applied urea or Super-Urea stabilized with N-butyl-thiophosphoric triamide and dicyandiamide urease and nitrification inhibitors (Spr urea or Spr super-U), and a no-manure or synthetic fertilizer application control (Control). Drystack manure and composted manure were applied to treatments at a target weight of 52 Mg ha^–1 and 33 Mg ha^–1 on a dry weight basis, respectively. All plots were disked following manure, compost, or synthetic fertilizer applications to incorporate treatments. Synthetic fertilizer application was based on pre-plant soil testing and University of Idaho guidelines for each crop [24]. As Fall compost treatments did not meet nutrient recommendations, synthetic fertilizer was supplemented accordingly. Additional information regarding site set-up is reported in Dungan et al. [26].

2.2. Calibration Data

Although the focus of this study was SOC stock, multiple parameters were used to assess and calibrate the DNDC model [21] for each site. For annual SOC stocks, soil samples were collected each fall before treatment application to a depth of 122 cm using a hydraulic soil probe (9100 Ag Probe, AMS Inc. American Falls, ID). Soil cores were separated into five segments (0–15, 15–30, 31–60, 61–91, and 92–122 cm) airdried and ground to pass a 2 mm sieve before SOC analysis by dichromate oxidation on a microplate spectrophotometer [23,27]. A correction factor of 1.33 was applied for incomplete oxidation of SOC [28]. For the purposes of this study, SOC contents were utilized from 0–30 cm as significant changes to SOC at lower depth intervals were not encountered over the duration of each study. Soil organic carbon contents were calculated on an area basis (Mg ha^–1) using estimated values of soil bulk density for 0–15 cm and 15–30 cm, of 1.28 g cm³ and 1.37 g cm³, respectively.

Crop yields were measured during mechanical harvest of research plots at each site; biomass-C calculated from dry combustion in a FlashEA1112 C/N analyzer is presented in lieu of dry matter mass (CE Elantech, Lakewood, NJ, USA). At CoverCrop, corn was harvested as silage and triticale was harvested as a forage. At GRACEnet, corn was harvested as silage, barley as grain with barley residue baled and removed, and alfalfa as a forage with 2–3 cuttings per year. At LT manure, wheat and barley were harvested as grain (with residue baled and removed); potato and sugarbeet were harvested as tubers/roots. Yield data were available up until 2019 at CoverCrop and LT manure and 2018 at GRACEnet.

In-season estimates of N mineralization were made using a buried bag method employed in the study region [29]. Briefly, soil was sampled with a bucket auger from the 0–30 and 30–60 cm depths within each plot, composited and then packed into polyethylene tubes and reset in core holes to facilitate aerobic in-field incubation. Polyethylene bags were periodically removed from the field during the growing season and analyzed for inorganic N by flow injection colorimetry using Quickchem methods 12-107-06-2-A and 12-107-04-1-B [30,31]. For model assessment and calibration, cumulative N mineralization was utilized. Logistical requirements constrained N mineralization estimates at some sites and treatments. At CoverCrop, data were available for 2018 and 2019 seasons for CT-M-F, CT-M-T, CT-NM-F, and CT-NM-T while only the 2019 season was available for MT-NM-F and MT-NM-T treatments; no data existed for MT-M-F and MT-M-T treatments. At GRACEnet, data were available for 2018 and 2019 seasons except for Fall compost and Spr super-U treatments where no data existed. At the LT manure site, data were available for all treatments and seasons.

Observed values of volumetric soil water content were obtained by time domain reflectometry (GRACEnet and LT Manure) or neutron probe (CoverCrop) at a depth of 0–15 cm. Soil water content data were available for all seasons and treatments at GRACEnet and CoverCrop. At LT manure 2013 to 2017 records were available excluding 18B, 36A, and 52B treatments which had records for 2015 to 2017. Lastly, AET estimates from ETIdaho, a multi-crop ET estimator for Idaho, were compared with DNDC AET estimates between planting and harvest dates for model calibration [32,33]. When a perennial was grown year-round, the annual AET was used. Although the ETIdaho estimated AET is most similar to DNDC estimated AET (by accounting for crop stresses and more complex shallow soil water dynamics), the ETIdaho database was last updated in 2016 which limited comparisons. As a result, where ETIdaho data were lacking, a more generalized crop specific AET estimate (Agrimet) was compared with DNDC AET. Calibration was performed solely on ETIdaho estimates; however, figures and statistics were updated to include the Agrimet data (2017 to 2019) to enhance comparison. Generally, the Agrimet AET estimates are slightly larger than ETIdaho estimates, however trends in water use are similar.

2.3. DNDC

The DNDC models began 1 year prior to the beginning of each study to permit fall treatment applications to occur. The DNDC model was initialized after selecting weather records and designating site management practices while allowing soil and crop parameters to be populated by default values where applicable. For example, selecting a soil texture results in default values of porosity, hydraulic conductivity, field capacity, and wilting point. Weather data were gathered from an Agrimet weather station within close proximity of research plots in Kimberly, Idaho (www.usbr.gov/pn/agrimet accessed on 10 January 2021). Soil parameters supplied under the default scenario were texture class, soil pH, and initial SOC content as obtained at 0–15 cm the fall each study began. Cropping parameters supplied included rotation information, fraction of residue left in field after a harvest event, planting and harvest dates, tillage events, fertilization conditions, manure or compost applications, and irrigation schedules.

DNDC model performance was evaluated using goodness of fit indicators: percent bias (PBIAS), Nash-Sutcliffe model efficiency coefficient (NSE), and model mean absolute error (MAE). Model PBIAS determines the tendency of simulated values to be larger or smaller than corresponding observed values; the optimal PBIAS value is 0 with positive and negative values indicating overestimation or underestimation bias, respectively. The NSE describes the predictive accuracy of a model ranging from negative infinity to 1, a perfect model; when NSE is equal to 0 the model has the same predictive ability as the mean of observations. Typically, NSE is a goodness of fit metric used in hydrologic simulations of streamflow where NSE ≥ 0.75 is considered good and NSE ≥ 0.36 < 0.75 is considered satisfactory agreement [34]. Model MAE is the calculated average of absolute errors between simulated and observed values.

$$PBIAS=100\frac{{{\displaystyle \sum }}_{i=1}^n\left(si{m}_i-ob{s}_i\right)}{{{\displaystyle \sum }}_{i=1}^n\left(ob{s}_i\right)} ; NSE=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(ob{s}_i-si{m}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2} ; MAE=\frac{{{\displaystyle \sum }}_{i=1}^n\left|si{m}_i-ob{s}_i\right|}{n}$$

To initiate calibration, parameter values were replaced with observed values where measurements or estimates existed for crop parameters (max biomass production and partitioned biomass C:N) and soil parameters (initial N concentration, bulk density, field capacity, wilting point, clay fraction, and hydraulic conductivity). Preliminary calibration was attempted on a singular synthetic fertilizer treatment at each site; however, it became apparent that a calibration performed in this way would not accurately represent treatments where manure was applied. Therefore, two treatments were considered at each site for the calibration. At CoverCrop, CT-NM-T and CT-M-T treatments were considered; at GRACEnet, Spr urea and Fall manure treatments were considered; at LT manure, Fertilizer and 52A treatments were considered. The trial and error method was used to manually calibrate model output at each of the three locations to observed values of volumetric soil water contents, biomass C, SOC, cumulative N mineralization, and estimated AET in that order and reiterated until model performance was deemed adequate or further benefit to goodness of fit indicators was not obtained. The priority of parameter consideration was based upon available empirical data and previous recommendations [35]. Utilization of the open source parameter estimation software (PEST; pesthomepage.org) was not feasible due to the graphical user interface of the DNDC model. Upon completion of the calibration, remaining treatments were initialized in the model and goodness of fit indicators produced. Although quasi perfect model performance was not an anticipated outcome of this study, substantial consideration was given to the calibration which is the product of numerous individual model runs. To reasonably constrain the length of this report, complete summary statistics and visualizations of model fit are reported as Supplementary Material. Instead, summary statistics are reported for all 3 locations for each parameter used in calibration; PBIAS and NSE were recalculated based on all treatments.

2.4. Forecasting

The calibrated models were used to forecast SOC levels under a subset of treatments up until year 2050. All treatments were forecasted at CoverCrop; at GRACEnet, Spr urea, Fall and Spr manure were forecasted; at LT manure 18A, 36B, 52A and Fertilizer treatments were forecasted. Model event schedules were modified in order to loop cropping rotations for the length of the forecasted scenario. For CoverCrop, the continuous corn cropping system permitted all events (2015–2019) to be looped. For GRACEnet, the first complete rotation of corn-barley-alfalfa-alfalfa-alfalfa (2013–2017) was looped. For LT Manure, two cycles of the wheat–potato–barley–sugarbeet rotation (2013–2020) was looped after repeating the 2016 sugarbeet year in 2020. Climate files used for each respective looping period were also repeated to reduce variability in the forecasted scenario. Precipitation in climate files was modified to reflect a projected 7.5% increase in regional precipitation by 2050 [36]. Similarly, temperature was modified according to “high” or “low” emissions projections resulting in 2.64 and 1.94 °C higher temperatures on average by 2050 and output for both the “high” and “low” emissions future climate scenarios are reported. Atmospheric CO₂ was held constant during the scenario as the default value of 350 ppmv was not altered during the “default” and calibrated scenarios.

2.5. Statistical Analysis

To assess treatment effects on SOC stocks, observed values of SOC were analyzed by ANOVA performed by site and year in R (version 3.6.3; [37]). ANOVA data were checked for normality and homogeneity of variance assumptions and transformed using the bestNormalize package in R as necessary [38]. Treatment means were separated using Tukey’s multiple range test using the agricolae package in R [39]. Non-transformed data are shown in tables and figures. To calculate the average annual rate of change in SOC, simple linear regression was used to plot SOC level as a function of year among treatments at each site, slope estimates and their confidence intervals were extracted from the models. At each site, further testing was conducted based on points of interest.

At CoverCrop, to assess whether any effect of tillage practice or winter triticale could be separated from the effect of manure application a linear mixed effects model fitting fixed factors of manure (M & NM), cover (F & T), tillage (CT & MT), and year (2015:2019), their interactions, and random intercept effects attributed to block, block:cover:manure, and block:tillage:cover:manure was constructed using the lme4 package in R [40]. Data were log-transformed to satisfy assumptions of homogeneity of variance and normality of residuals before subsequent ANOVA and contrast testing was conducted to explain significant effects.

At GRACEnet, data were subset to manure applied treatments and a linear mixed effects model fitting timing (Fall & Spring), year (2012:2019), their interaction, and a random intercept effect of block was constructed to determine whether manure application timing (Fall & Spring) affected SOC accumulation. To determine if including 3 years of alfalfa in the rotation affected SOC stocks, data were subset to years before and directly proceeding the 3 years of alfalfa growth (2014 & 2017) before a Welch’s two sample paired t-test was conducted to assess if the difference in SOC means in 2014 and 2017 was different from 0.

At LT manure, a linear mixed effects model was fit with application rate (18, 36, 52 Mg ha^–1), frequency (annual, biennial), year (2012:2018), their interactions, and a random intercept effect of block to assess the difference, if any, in SOC accumulation between annual and biennial manure applications. Data were transformed to satisfy homogeneity of variance and residual normality assumptions before an ANOVA was conducted to assess the significance of model parameters. Here, contrast testing was conducted to assess the difference, if any, in SOC accumulation of 18A and 36B treatments.

The DNDC models were assessed in a “default” and calibrated capacity using goodness of fit indicators discussed above. Overall location (CoverCrop, GRACEnet, LT manure) statistics were recalculated using data from all simulated treatments, in the case of MAE this is the mean of all treatments. Where applicable all tests were considered significant at the 0.05 level; in one explicit instance, a marginally significant effect was considered at the 0.10 level.

3. Results and Discussion

3.1. Annual SOC Stocks

Soil organic carbon contents of the surface soil (0–30 cm) was measured annually at each of three locations from the beginning of each study to 2019. There was no known history of manure application before the commencement of each study. As a result, a wide range in SOC contents was observed between all three sites from 25.1 Mg kg^–1 to 89.6 Mg kg^–1 (

Table 1. Soil organic carbon contents obtained from fall soil sampling from 0–30 cm. Soil organic carbon contents were calculated on an area basis using assumed bulk densities of 1.28 g cm³ and 1.37 g cm³ for 0–15 cm and 15–30 cm depths, respectively ^1,2.

	Total Soil Organic Carbon 0–30 cm
Treatment	2012	2013	2014	2015	2016	2017	2018	2019
	––––––––––––––––––––––––––– Mg ha−1 –––––––––––––––––––––––––––
CoverCrop
CT-M-F	-	-	-	27.2 a	37.5 ab	49.3 a	57.8 a	66.5 a
CT-M-T	-	-	-	28.8 a	39.0 a	51.0 a	50.8 a	66.8 a
CT-NM-F	-	-	-	26.2 a	28.8 c	29.8 b	29.0 b	35.3 b
CT-NM-T	-	-	-	26.7 a	28.6 c	31.4 b	29.6 b	35.5 b
MT-M-F	-	-	-	26.8 a	37.3 ab	47.2 a	59.8 a	69.9 a
MT-M-T	-	-	-	25.2 a	37.1 ab	48.9 a	58.5 a	89.6 a
MT-NM-F	-	-	-	25.9 a	30.8 bc	29.2 b	29.2 b	33.0 b
MT-NM-T	-	-	-	25.1 a	29.2 c	27.9 b	30.2 b	34.4 b
p				0.63	<0.001	<0.001	<0.001	<0.001
GRACEnet
Control	38.5 a	38.1 bc	37.9 c	42.5 b	40.5 b	39.6 b	41.8 b	42.2 c
Fall compost	39.2 a	42.0 abc	44.1 b	43.8 b	43.8 ab	42.8 b	46.9 ab	49.0 bc
Fall manure	40.5 a	47.5 ab	54.2 a	60.7 a	52.5 a	49.5 a	53.8 a	72.2 a
Spr manure	39.1 a	48.7 a	52.1 a	60.6 a	55.2 a	52.4 a	54.0 a	61.1 ab
Spr super-U	44.3 a	36.6 c	40.3 bc	41.9 b	37.3 b	40.2 b	42.9 b	41.9 c
Spr urea	37.8 a	41.4 abc	39.0 c	40.3 b	39.1 b	40.9 b	41.9 b	43.0 c
p	0.08	0.01	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001
LT manure
18A	40.8 a	36.6 bcd	42.5 abc	51.5 abc	38.9 bc	45.3 bcd	41.6 b	55.4 bc
18B	41.4 a	39.4 abcd	36.8 bc	42.4 cd	37.5 bc	38.6 de	37.8 b	48.8 c
36A	43.2 a	40.3 abcd	45.9 ab	58.3 ab	54.8 a	57.7 ab	45.5 ab	67.5 ab
36B	40.4 a	43.9 abc	40.9 bc	45.3 bcd	38.4 bc	44.2 cd	41.2 b	51.9 c
52A	38.7 a	49.9 a	53.9 a	69.7 a	57.5 a	67.9 a	59.8 a	83.7 a
52B	40.9 a	44.9 ab	40.0 bc	69.2 a	47.0 ab	54.0 abc	43.1 ab	67.9 ab
Control	40.9 a	33.5 d	35.3 c	37.4 d	33.5 c	34.0 e	37.7 b	39.4 d
Fertilizer	39.3 a	35.3 cd	35.4 c	37.8 d	34.7 c	34.5 e	37.9 b	41.5 d
p	0.29	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001	<0.001

¹ Values within each column and site containing the same letter are not significantly different at the α = 0.05 level by Tukey’s multiple range test. The p-values at the bottom of each column reflect the F-statistic of the overall ANOVA. ² Cover Crop treatments: CT-M-F, conventional tillage-manure application-fallow; CT-M-T, conventional tillage-manure application-triticale; CT-NM-F, conventional tillage-no manure-fallow; CT-NM-T, conventional tillage-no manure-triticale; MT-M-F, Minimum-till-manure application-fallow; MT-M-T, Minimum-till-manure application-triticale; MT-NM-F, Minimum-till-no manure-fallow; MT-NM-T, Minimum-till-no manure-triticale. Treatments with manure were applied annually at a target rate of 52 Mg ha⁻¹ on a dry weight basis. GRACEnet treatments: Control, no synthetic fertilizer or manure; Fall compost, composted dairy manure applied in the Fall; Fall manure, dairy manure applied in the Fall; Spr manure, dairy manure applied in the Spring; Spr super-U, Super-U applied annually in the spring based on soil test N; Urea, Urea applied annually in the spring based on soil test N. Compost and manure applications were made on a dry weight basis according to crop rotation at target application rates of 33 Mg ha^–1 and 52 Mg ha^–1, respectively. LT manure treatments: 18A, dairy manure applied annually at a rate of 18 Mg ha^–1; 18B, dairy manure applied biennially at a rate of 18 Mg ha^–1; 36A, dairy manure applied annually at a rate of 36 Mg ha^–1; 36B, dairy manure applied biennially at a rate of 36 Mg ha^–1; 52A, dairy manure applied annually at a rate of 52 Mg ha^–1; 52B, dairy manure applied biennially at a rate of 52 Mg ha^–1; Control, no synthetic fertilizer or manure; Fertilizer, synthetic fertilizer. All target manure application rates are on a dry weight basis.

).

To assist interpretation, using assumed bulk densities of 1.28 g cm³ and 1.37 g cm³ for 0–15 cm and 15–30 cm depths, respectively, 20 Mg ha^–1 is slightly less than 0.5% SOC. Average SOC content at the CoverCrop site was numerically lower when the study was initiated in 2015 (26.5 Mg ha^–1) than when either GRACEnet (39.9 Mg ha^–1) or LT manure (40.7 Mg ha^–1) studies were initiated in 2012. This can be attributed to the CoverCrop study being located on the head of a historically furrow irrigated field, as these locations are highly erosive [41,42]. At the CoverCrop site, by 2017 SOC content was universally significantly higher where manure was applied, and this persisted through 2019. Conversely, where manure had not been applied SOC increased by approximately 30% by 2019; however, any effect of triticale cover or tillage practice was not apparent. At the GRACEnet location, in no year were SOC stocks different between synthetic fertilizer treatments and control treatments. By 2019, SOC of the fall manure treatment was ~18% higher than, but not significantly different from the spring applied manure treatment. Curiously, any significant difference between fall and spring manure applications was not apparent from fall 2013 to 2015 even though more manure-C had been applied under the Fall manure treatment during this time (

Table 2. Manure carbon inputs and biomass carbon removed from yield and residues where applicable ¹.

	Applied Manure-C									Biomass-C Removed
	2012	2013	2014	2015	2016	2017	2018	2019	Cumulative	2012	2013	2014	2015	2016	2017	2018	2019	Cumulative
CoverCrop	––––––––––––––––––––––––––– Mg ha–1 –––––––––––––––––––––––––––
CT-M-F	-	-	-	14.92	16.70	9.17	8.88	10.84	60.51	-	-	-	-	8.26	9.12	7.46	7.19	32.03
CT-M-T	-	-	-	14.30	17.16	9.25	8.19	10.83	59.73	-	-	-	-	8.06	13.35	8.03	9.61	39.05
CT-NM-F	-	-	-	-	-	-	-	-	-	-	-	-	-	8.64	8.97	7.22	6.86	31.69
CT-NM-T	-	-	-	-	-	-	-	-	-	-	-	-	-	8.04	9.77	7.80	7.36	32.98
MT-M-F	-	-	-	14.72	16.94	9.00	8.18	10.85	59.68	-	-	-	-	8.47	8.96	7.64	7.26	32.32
MT-M-T	-	-	-	15.50	17.45	9.89	8.56	10.84	62.24	-	-	-	-	8.34	12.54	8.24	9.45	38.57
MT-NM-F	-	-	-	-	-	-	-	-	-	-	-	-	-	8.06	7.87	7.65	6.76	30.33
MT-NM-T	-	-	-	-	-	-	-	-	-	-	-	-	-	8.65	8.67	7.61	6.98	31.91
GRACEnet
Control	-	-	-	-	-	-	-	-	-	-	9.29	3.72	2.79	5.67	6.08	10.84	-	38.39
Fall compost	2.98	3.38	-	-	-	3.30	1.78	-	11.45	-	8.28	4.29	3.30	6.28	6.88	9.13	-	38.16
Fall manure	15.97	15.45	-	-	-	9.26	7.99	-	48.67	-	8.13	6.21	3.41	6.34	7.12	10.54	-	41.75
Spr manure	-	7.93	13.22	-	-	-	5.10	10.06	36.30	-	8.18	4.99	3.56	6.20	6.78	8.91	-	38.61
Spr super-U	-	-	-	-	-	-	-	-	-	-	8.91	5.77	3.24	6.24	6.54	10.99	-	41.69
Spr urea	-	-	-	-	-	-	-	-	-	-	9.00	6.06	3.17	6.08	6.80	11.69	-	42.78
LT manure
18A	4.94	6.40	3.75	4.02	3.37	2.80	2.17	2.10	29.54	-	2.90	4.01	6.20	7.26	5.64	4.18	6.63	36.81
18B	4.94	-	3.75	-	3.37	-	2.17	-	14.22	-	2.66	4.02	6.14	7.96	5.58	3.80	6.29	36.46
36A	9.94	12.89	7.54	8.10	6.78	5.63	4.36	4.23	59.46	-	3.01	3.78	6.38	7.27	5.13	3.62	6.53	35.71
36B	9.94	-	7.54	-	6.78	-	4.36	-	28.61	-	2.55	4.00	6.24	7.76	5.82	4.34	6.70	37.41
52A	14.87	19.30	11.29	12.12	10.14	8.43	6.53	6.33	89.00	-	2.83	3.26	6.92	8.20	5.57	3.64	5.61	36.03
52B	14.87	-	11.29	-	10.14	-	6.53	-	42.83	-	2.52	4.30	6.79	6.96	5.83	4.59	6.05	37.04
Control	-	-	-	-	-	-	-	-	-	-	2.40	3.61	3.90	6.63	2.78	3.47	2.95	25.74
Fertilizer	-	-	-	-	-	-	-	-	-	-	2.43	4.29	5.65	8.08	5.82	3.89	5.72	35.88

¹Cover Crop treatments: CT-M-F, conventional tillage-manure application-fallow; CT-M-T, conventional tillage-manure application-triticale; CT-NM-F, conventional tillage-no manure-fallow; CT-NM-T, conventional tillage-no manure-triticale; MT-M-F, Minimum-till-manure application-fallow; MT-M-T, Minimum-till-manure application-triticale; MT-NM-F, Minimum-till-no manure-fallow; MT-NM-T, Minimum-till-no manure-triticale. Treatments with manure were applied annually at a target rate of 52 Mg ha⁻¹ on a dry weight basis. GRACEnet treatments: Control, no synthetic fertilizer or manure; Fall compost, composted dairy manure applied in the Fall; Fall manure, dairy manure applied in the Fall; Spr manure, dairy manure applied in the Spring; Spr super-U, Super-U applied annually in the spring based on soil test N; Urea, Urea applied annually in the spring based on soil test N. Compost and manure applications were made on a dry weight basis according to crop rotation at target application rates of 33 Mg ha^–1 and 52 Mg ha^–1, respectively. LT manure treatments: 18A, dairy manure applied annually at a rate of 18 Mg ha^–1; 18B, dairy manure applied biennially at a rate of 18 Mg ha^–1; 36A, dairy manure applied annually at a rate of 36 Mg ha^–1; 36B, dairy manure applied biennially at a rate of 36 Mg ha^–1; 52A, dairy manure applied annually at a rate of 52 Mg ha^–1; 52B, dairy manure applied biennially at a rate of 52 Mg ha^–1; Control, no synthetic fertilizer or manure; Fertilizer, synthetic fertilizer. All target manure application rates are on a dry weight basis.

).

As perhaps expected, the compost application was numerically between the manure applications and the remaining treatments, but not significantly different from either. When observing SOC stocks of the LT manure site, one trend encountered was an apparent reduction in SOC values the fall succeeding potato and sugarbeet harvest. In 2016 and 2018, mean SOC levels decreased an average of 15.7% and 6.2% relative to the respective preceding year. This could be a product of homogenization of soil depths during tuber/root harvest; alternatively, low residue C input, increased decomposition due to low residue C:N, and removal of belowground biomass during these crops has been attributed to decreases in SOC [43,44,45]. The current study anecdotally supports a degradative coefficient [46] for potato and sugar beet; though, our data is limited. Nonetheless, when tracking SOC content over time it could be advantageous to standardize the point of comparison in a crop rotation that includes tuber/root crops. In no year was there a significant difference in SOC content between the Control and Fertilizer treatments (Table 1). As anticipated, by 2019 the annual manure applications and increasing application rates resulted in numerically higher SOC levels than biennial and lower application rate counterparts, nevertheless the differences were not always significant.

3.2. Average Rate of Change

Average rate of SOC change in each treatment was determined using simple linear regression at each site (

Table 3. Slope estimate and confidence interval for simple linear regressions fit to observed annual mean soil organic carbon contents, 0–30 cm ¹.

	Slope
Treatment	Estimate	C.I. 2.5%	C.I. 97.5%	p > F
CoverCrop	––––––––––––––––––––––– Mg ha–1 year−1 ––––––––––––––––––––––
CT-M-F	9.9	8.1	11.7	<0.001
CT-M-T	8.8	7.0	10.6	<0.001
CT-NM-F	1.8	0.8	2.9	<0.010
CT-NM-T	1.9	0.8	2.9	<0.010
MT-M-F	10.9	8.3	13.4	<0.001
MT-M-T	15.0	12.2	17.8	<0.001
MT-NM-F	1.2	0.3	2.2	<0.050
MT-NM-T	1.9	1.2	2.7	<0.001
GRACEnet
Control	0.6	<0.1	1.1	<0.050
Fall compost	1.1	0.6	1.6	<0.001
Fall manure	2.8	1.4	4.1	<0.001
Spr manure	2.1	1.3	2.9	<0.001
Spr super-U	0.1	–0.5	0.8	0.719
Spr urea	0.5	0.1	0.9	<0.050
LT manure
18A	1.5	0.5	2.4	<0.010
18B	0.5	–0.3	1.3	0.188
36A	2.7	1.3	4.1	<0.001
36B	0.8	0.0	1.7	0.061
52A	4.7	3.1	6.3	<0.001
52B	2.4	0.7	4.0	<0.010
Control	0.0	–0.6	0.7	0.905
Fertilizer	0.3	–0.2	0.7	0.273

¹Cover Crop treatments: CT-M-F, conventional tillage-manure application-fallow; CT-M-T, conventional tillage-manure application-triticale; CT-NM-F, conventional tillage-no manure-fallow; CT-NM-T, conventional tillage-no manure-triticale; MT-M-F, Minimum-till-manure application-fallow; MT-M-T, Minimum-till-manure application-triticale; MT-NM-F, Minimum-till-no manure-fallow; MT-NM-T, Minimum-till-no manure-triticale. Treatments with manure were applied annually at a target rate of 52 Mg ha^–1 on a dry weight basis. GRACEnet treatments: Control, no synthetic fertilizer or manure; Fall compost, composted dairy manure applied in the Fall; Fall manure, dairy manure applied in the Fall; Spr manure, dairy manure applied in the Spring; Spr super-U, Super-U applied annually in the spring based on soil test N; Urea, Urea applied annually in the spring based on soil test N. Compost and manure applications were made on a dry weight basis according to crop rotation at target application rates of 33 Mg ha^–1 and 52 Mg ha^–1, respectively. LT manure treatments: 18A, dairy manure applied annually at a rate of 18 Mg ha^–1; 18B, dairy manure applied biennially at a rate of 18 Mg ha^–1; 36A, dairy manure applied annually at a rate of 36 Mg ha^–1; 36B, dairy manure applied biennially at a rate of 36 Mg ha^–1; 52A, dairy manure applied annually at a rate of 52 Mg ha^–1; 52B, dairy manure applied biennially at a rate of 52 Mg ha^–1; Control, no synthetic fertilizer or manure; Fertilizer, synthetic fertilizer. All target manure application rates are on a dry weight basis.

).

At the CoverCrop location, slopes of the regressions were exceptionally high where manure was applied, indicating annual SOC accumulations of ~ 10 Mg ha^–1 year^–1. This is considerably higher than studies of comparative lengths of manure application in the literature [47,48,49] and at the other locations of the present study. However, considering the relatively extreme manure application rate at CoverCrop (52 Mg ha^–1 year^–1), the cumulative manure-C input amassed quickly (Table 2). Additional SOC accumulation may be related to the shift to a continuous corn cropping system thought to be beneficial in some cases due to an increase in residue-C input [50,51,52]. Corn stover retention is estimated to increase SOC stocks at a rate of 0.41 Mg ha^–1 year^–1 relative to baseline values [53] but was removed for silage in the present study. Elsewhere, Bolinder et al. [54] attributed a ~0.3 Mg ha^–1 year^–1 SOC increase mainly to corn root derived C. Therefore, average annual C additions of manure applied plots at CoverCrop may approach 12 Mg ha^–1 year^–1. If we assume a moderate estimated SOC maintenance C addition rate of 1.9 Mg ha^–1 year^–1 [55], our estimate of SOC accumulation is conceivable. Additional explanation may be attributed to an additive effect of a less erosive irrigation regime (furrow to sprinkler) which was adopted when this study was initiated. Accumulation of SOC where no manure was applied was far lower, nonetheless considerable given the SOC accumulations at GRACEnet and LT manure locations (Table 3). As a result, this study indicates that SOC can be accrued rapidly from a diminished point by high manure-C input, although there are other considerations for high manure application rates in this region such as soil P and soluble salt accumulation [56].

At GRACEnet, the slope of the regression for the Spr super-U treatment was not significantly different from 0; consequently, SOC had not meaningfully changed since initial 2012 levels (Table 3). We propose that accumulation rates were lower at GRACEnet relative to CoverCrop for two reasons. First, CoverCrop received ~40% more manure-C over the span of 3 fewer years (Table 2). Second, initial SOC level was considerably higher at GRACEnet and responses to C inputs are related to initial SOC levels and effective equilibriums [57,58].

The slope estimates at LT manure were not significantly different from 0 in Control, Fertilizer, 18B and 36B treatments; suggesting in this cropping system manure-C input of at least 18 Mg ha^–1 annually or 36 Mg ha^–1 biennially is necessary for accrual of SOC (Table 3). As the GRACEnet and LT manure studies were initiated in 2012 at approximately equivalent SOC levels on very similar soils, they are relevant to compare. The control and synthetic fertilizer (Spr urea for GRACEnet) treatments appear to be accumulating SOC in the GRACEnet study but not at LT manure (Table 3). This was attributed to both a lesser amount of soil disturbance and potentially higher residue C input under the GRACEnet cropping rotation; residue C under potato or sugarbeet has been estimated as 2 to 3 times lower than residue C from cereal crops [59]. Our estimates of SOC accumulation can be compared with others but is limited by the scarcity of SOC records in potato or sugarbeet rotations receiving manure. Our own estimate using data from Moulin et al. [60] suggests SOC accumulation of 2.7 Mg ha^–1 year^–1 under a composted cattle manure application rate of ~ 15 Mg ha^–1 year^–1 in a 5-year bean-potato succession. Elsewhere, data from Miao et al. [61] indicated SOC accrual (~ 3 Mg ha^–1 year^–1) under a similar application rate (14 Mg ha^–1 year^–1) in a corn soybean succession; both studies estimates were higher than in the 18A treatment of the present study (1.5 Mg ha^–1 year^–1, Table 3).

3.3. Specific Management Questions

At CoverCrop, a linear mixed effects model was used to determine any significant effect of tillage practice or winter cover (triticale) on SOC stocks irrespective of manure application. An ANOVA of the model indicated significant main effects of manure application and year (both p < 0.0001), no effect of winter cover (p = 0.51) and an interaction effect between manure, tillage, and year (p < 0.05). Contrast testing indicated that tillage practice effected SOC levels only in 2019, solely where manure was applied (p < 0.01). This appears to be the result of the large individual mean measured under the MT-M-T treatment in 2019 relative to other manure applied treatments (Table 1). Consequently, we do not believe there is conclusive evidence that tillage practice or winter triticale cover has affected SOC accumulation in this system. Meta-analyses in the literature are often contradictory on the effect of tillage intensity on total SOC stock. Some are in agreement with the present study [62] others are not [63,64], but most report differences in SOC depth distribution and utilize only studies beyond a selected duration length. The CoverCrop study duration is ~ 5 years and would have been included in each of the mentioned meta-analyses. Conversely, meta-analyses considering the impact of cover crops on SOC stock consistently report accumulation at rates of 0.2 to 0.6 Mg ha^–1 year^–1 [65,66,67,68]. The present study encountered limited triticale growth where manure had not been applied before corn had to be planted (Figure S6); additionally, triticale was harvested as forage instead of returning to the soil as a green manure. Both are possible explanations for why triticale cover did not have a significant effect on SOC stock.

A linear mixed effects model was used to determine the difference, if any, in SOC accumulation between fall and spring applied manure at the GRACEnet location. Results indicated a significant main effect of year (p < 0.0001), no effect of timing (p = 0.41), and a timing (fall vs. spring) by year interaction that could be considered marginally significant (p < 0.1). Contrast testing of this interaction indicated the only difference between the effect of timing on SOC levels was in 2019 (p < 0.001). Considering individual years, the SOC stock of Fall manure and Spr manure treatments differed <6% from 2012–2018. Therefore, we find that there is insufficient evidence to state that there is a significant difference in SOC accumulation between fall and spring applications of manure. There are a limited number of comparative studies of manure application timing on SOC stock. One study reported higher SOC content in 1 of 2 years under fall application relative to spring application [69]; Ahmed et al. [70] reported higher organic matter under injected spring applications relative to fall applications after 6 years. Theoretically, if total manure-C input is equal it is likely there is not a substantial difference between timing as SOC stock changes are related to cumulative total C inputs [71]. One reason this may not be true is if differences in nitrogen release and retention between timings results in increases in above and below ground crop biomass. In this respect, there may be a slight preference for spring applications as fall application leaves an elongated window for nitrogen losses, although, precipitation dynamics must also be considered [70,72,73]. In the region of study, total precipitation is low (~24 cm year^–1) and characterized by low intensity rainfall events [74]; that in addition to well managed sprinkler irrigation systems suggests minimal leaching of N.

In the present study, the effect of including alfalfa in the GRACEnet dairy forage rotation on SOC stock was also of interest. A paired Welch’s t-test was used to determine any difference in SOC stocks before and after alfalfa was grown in the rotation (2014 and 2017). The estimated stock difference (0.18 ± 0.68 Mg ha^–1) was not significantly different from 0 (p = 0.59). Nonetheless, as the Control treatment accumulated SOC from 2012 to 2019 at GRACEnet but not at LT manure, the overall effect of the dairy forage rotation on SOC can be considered positive (Table 3). In comparable studies, alfalfa has been credited for improving SOC when included in rotations due to a reduction in tillage frequency and high residue C input [75,76].

At the LT manure study, a linear mixed effects model was used to determine the difference, if any, between manure application frequency (annual vs. biennial). Model fixed effects showed no significant difference (p = 0.75) between application frequency on SOC stock. Another point of comparison was made after an ANOVA of the model indicated a significant (p < 0.05) interaction between application rate, frequency, and year. Contrast testing of frequency and rate combinations receiving equivalent quantities of manure (18A and 36B) suggested differences in SOC stock only in the first year of manure application, thus annual application at rate 1X is approximately equivalent to biennial application at rate 2X. These findings are in agreement with Eghball [77] who reported no difference in SOC stock after 4 years of annual and biennial application of cattle manure at rate 1X and rate 2X, respectively.

3.4. DNDC Modeling

Overall model performance was dependent upon the parameter of interest and study location. The overall model goodness of fit statistics for both the default and calibrated models are reported in

Table 4. Overall denitrification decomposition (DNDC) model goodness of fit statistics for each location. Overall location statistics were calculated using data from all treatments, in the case of mean absolute error (MAE), this is the mean of all treatments ¹.

	PBIAS		NSE		MAE
	Default	Calibrated	Default	Calibrated	Default	Calibrated
	––––––––––––––––– Soil Organic Carbon, Mg ha−1 –––––––––––––––––
CoverCrop	−9.8	−14.4	0.59	0.48	5.08	6.06
GRACEnet	−6.9	−3.9	0.45	0.76	4.74	2.75
LT manure	0.6	2.9	0.74	0.73	4.12	4.36
	–––––––––––––––– Biomass Carbon (yield), Mg ha−1 ––––––––––––––––
CoverCrop	4.1	5.0	0.84	0.92	0.99	0.58
GRACEnet	−32.6	8.3	−0.69	0.70	2.62	1.13
LT manure	−14.4	−4.4	0.73	0.85	0.67	0.49
	–––––––––––––––– Nitrogen Mineralization, kg ha−1 ––––––––––––––––
CoverCrop	−10.0	20.2	0.37	0.13	51	61
GRACEnet	−42.4	−33.7	0.33	0.54	59	48
LT manure	−30.9	6.2	0.70	0.92	121	59
	––––––––––––––––––––––– Soil Water, Θv –––––––––––––––––––––––
CoverCrop	34.1	27.9	−4.72	−3.49	0.09	0.08
GRACEnet	46.1	9.4	−2.03	−0.03	0.11	0.06
LT manure	14.3	4.2	−0.73	−1.07	0.08	0.09
	–––––––––––– Estimated Actual Evapotranspiration, mm ––––––––––––
CoverCrop	−44.5	−43.7	−1.17	−1.11	244	240
GRACEnet	−37.4	−34.0	−1.92	−1.55	303	275
LT manure	−14.3	−7.9	0.25	0.59	150	80

¹Cover Crop treatments: CT-M-F, conventional tillage-manure application-fallow; CT-M-T, conventional tillage-manure application-triticale; CT-NM-F, conventional tillage-no manure-fallow; CT-NM-T, conventional tillage-no manure-triticale; MT-M-F, Minimum-till-manure application-fallow; MT-M-T, Minimum-till-manure application-triticale; MT-NM-F, Minimum-till-no manure-fallow; MT-NM-T, Minimum-till-no manure-triticale. Treatments with manure were applied annually at a target rate of 52 Mg ha^–1 on a dry weight basis. GRACEnet treatments: Control, no synthetic fertilizer or manure; Fall compost, composted dairy manure applied in the Fall; Fall manure, dairy manure applied in the Fall; Spr manure, dairy manure applied in the Spring; Spr super-U, Super-U applied annually in the spring based on soil test N; Urea, Urea applied annually in the spring based on soil test N. Compost and manure applications were made on a dry weight basis according to crop rotation at target application rates of 33 Mg ha^–1 and 52 Mg ha^–1, respectively. LT manure treatments: 18A, dairy manure applied annually at a rate of 18 Mg ha^–1; 18B, dairy manure applied biennially at a rate of 18 Mg ha^–1; 36A, dairy manure applied annually at a rate of 36 Mg ha^–1; 36B, dairy manure applied biennially at a rate of 36 Mg ha^–1; 52A, dairy manure applied annually at a rate of 52 Mg ha^–1; 52B, dairy manure applied biennially at a rate of 52 Mg ha^–1; Control, no synthetic fertilizer or manure; Fertilizer, synthetic fertilizer. All target manure application rates are on a dry weight basis.

.

Generally, SOC values were simulated well by the default DNDC model. Under the default capacity, model PBIAS was ± 10% and NSE ranged from 0.45 to 0.74 (Table 4). The MAE indicated on average, SOC estimates were within 4.65 Mg ha^–1, approximately 1.2 mg kg^–1, of their measured values. Calibration attempts resulted in a slight deterioration of SOC fit statistics for CoverCrop and LT manure locations and a slight improvement at the GRACEnet location (Table 4). Unsurprisingly, individual treatment MAE indicated that both the default and calibrated DNDC models encountered the largest absolute error where manure was applied (Table S2). However, a tendency for underprediction of SOC when manure was applied was only encountered at the CoverCrop location which had unusually high SOC accumulation (Table 1 and Figure 1).

The MAE was considerably higher under minimum-till treatments (Table S2); subsequent examination of output files indicated that decomposition of C was questionably high (≥100 ≤ 500 kg^–1 ha^–1 day^–1) when manure was unincorporated. We suggest that DNDC modeled decomposition of manure applied C was excessively rapid when not incorporated (Figure 1). This presents an opportunity for model progression in future iterations; in a report on DNDC model development, Li et al. [78]:

“…a change in manure application depth can simultaneously alter the soil temperature, moisture, pH, Eh, and concentrations of DOC, NH₄⁺ or NO₃⁻. These changes will simultaneously and collectively affect the rates of decomposition, nitrification and denitrification occurring in the manure-amended soil that eventually alters the emissions of CO₂, N₂O and NH₃ from the soil.”

This result is also demonstrative of the utility of pairing empirical data with modeled output, even though adequate simulation of SOC by DNDC is typically reported [79,80]. Overall, our results indicate DNDC simulated SOC values adequately over the duration of our research studies in a default capacity considering the diverse management practices simulated, although calibration can improve model performance.

Model simulation of crop yield was assessed through biomass-C. The default model performance was highly dependent on the cropping system being simulated. The continuous corn system of CoverCrop had a PBIAS of 4.1% in default capacity and 5.0% after calibration, the calibration markedly improved simulation of triticale growth improving NSE from 0.84 to 0.92 and MAE from ~ 1 Mg C ha^–1 to 0.58 Mg C ha^–1 (Table 4 and Figure 2).

The commercial rotation of LT manure (wheat–potato–barley–sugarbeet) was simulated acceptably in the default capacity, PBIAS= –14.4%, NSE = 0.73, MAE = 0.67 Mg C ha^–1, but improved under calibration, mostly by increasing sugarbeet yield where manure was not applied (Table 4, Figures S26 and S27). The commercial forage rotation of GRACEnet poorly simulated crop yields under the default condition, PBIAS = –32.6%, NSE = –0.69, MAE = 2.62 Mg C ha^–1, as result of consistent underprediction of alfalfa yield and a failure of scheduled cuttings to occur in 2017 (Table 4 and Figure S16). An explanation for the failed cuttings in 2017 was not ascertained, however all cuttings occurred after calibration (Figure S17). Calibration of the GRACEnet location improved upon the default model, PBIAS = 8.3%, NSE = 0.70, MAE = 1.13 Mg C ha^–1, mostly by changing alfalfa growth parameters similar to those reported by He et al. [81]. Under the calibration, alfalfa growth simulation was overestimated in 2015 (Figure S17). A flaw forcing perennial crops to accumulate thermal degree days starting on 1 January of each year results in overestimation of yields when simulating spring planted perennials; subsequent model iterations should permit an exception to this rule. An additional peculiarity encountered under the default condition was an apparent lack of restriction on alfalfa root growth that resulted in an abrupt increase in SOC when alfalfa was terminated and root biomass was transferred to SOC pools (Figure S18). Alteration of the alfalfa growth parameters in the calibration addressed this concern (Figure S19). He et al. [81], reported “fair” simulation (NSE > 0 and nRSME > 30%) of alfalfa yield before improving simulation of winterkill effects, but did not report problems with root growth or spring planting. It is our view that biomass C was adequately simulated by the default model parameters if the alfalfa error is overlooked, yet attentiveness to model output is advised.

The DNDC model consistently underpredicted cumulative N mineralization in the default capacity where manure was not applied (Figure 3).

This was especially evident at the LT manure location as goodness of fit statistics improved with increasing application rates and were more favorable under annual as opposed to biennial applications (Figure 3, Table S3 and Figure S30). In response, decomposition rates were increased for the litter and humus pools under the calibration at all locations. The modeled N mineralization at the LT manure location was improved from the default scenario, PBIAS from –30.9% to 6.2%, NSE from 0.70 to 0.92, and MAE from 121 kg ha^–1 to 59 kg ha^–1 (Table 4 and Figure S31). Although decomposition rates were increased by the same margin at the GRACEnet location, fit was only modestly improved. At CoverCrop, overall goodness of fit indicators declined after calibration as only CT-NM-F and CT-NM-T treatments were improved (Table 4, Table S3). However, the quantity of comparable data was low at this location (Figures S10 and S11). Similar evaluations of the DNDC model have reported underestimation of N mineralization when contrasting various soil N pools [82,83], more so where N inputs were low [84,85]. Smith et al. [82] recommended revision of the mineralization module within DNDC if improvements to the hydraulic simulation could be completed beforehand. Krobel et al. [83] suggested that while increasing the decomposition rates would ameliorate underprediction of soil N pools, it may be detrimental to the soil C budget. In our calibration, decomposition rates were increased without degrading SOC simulation in 2 of 3 locations; although, changes to cropping parameters were required (Table 4).

The DNDC model was unable to acceptably simulate 0–15 cm soil water contents at all three locations under the default and initial calibration scenarios (Table 4 and Figure 4).

The calibration improved goodness of fit statistics at all locations, with reasonable PBIAS at GRACEnet (9.4%) and LT manure (4.2%); nevertheless, the NSE remained below 0 suggesting the mean soil water content was a better predictor of temporal soil water dynamics than the DNDC model (Table 4). A similar result was obtained by Jiang et al. [86] who reported reasonable simulation of soil water within an acceptable range (PIBAS ~ 3%) but worse simulated trends. Presently, soil water content was overestimated at all locations (PBIAS= 9–46%); elsewhere, comparable ranges of overestimation have been reported by Abdalla et al. [84] (13–30%) and Uzoma et al. [87] (13–26%). There is a consensus amongst researchers that a simplified approach to hydrologic simulation likely accounts for errors in modeled soil water content [88]. An additional consideration is that soil water content is a dynamic property, especially in irrigation driven agriculture, meaning it is likely difficult for models to correctly simulate values on specific days. The DNDC model uses a “tipping bucket”/cascade approach to model soil water movement; water entering the profile must saturate each soil layer before advancing to successive layers. While computationally less demanding, this method cannot account for movement of water in the vadose zone. A more precise approach is given by numerical solution to the Richards equation for movement of water in unsaturated media but has not been integrated in DNDC at this time. Additionally, water passing below a depth of 1 m is considered leached [89]. The present study indicates that the computationally inexpensive cascade approach is not suitable for modeling soil water trends in semi-arid irrigated croplands. Moreover, as water filled pore space is integral to other sub-models (nitrification, denitrification, decomposition), the accuracy of modelled emissions of GHG’s (CH₄, CO₂, N₂O) should be further evaluated for semi-arid regions.

Generally, the DNDC model estimated AET did not agree with ETIdaho estimates and was comparatively under predictive of water use (Figure 5).

The DNDC model estimated crop AET was simulated poorly at CoverCrop and GRACEnet locations under the default scenario with a PBIAS of –36.8% and –39.7%, respectively (Table 4). At CoverCrop, this appeared to be largely related to triticale water use affecting AET under the corn crop, despite soil water being overestimated by DNDC (Figure 5, Table S4 and Figure S14). One occurrence of near perfect agreement was seen in year 2015 at CoverCrop, probably as the field was fallowed in this year and no crop was simulated. At GRACEnet, estimated AET was farthest from ETIdaho estimates under years alfalfa was grown, likely exacerbating the overall NSE (–1.92) which was indicative of poor simulation of trends in water use by each crop (Table 4, Figures S24 and S25). Errors in alfalfa biomass simulation likely contributed to poor AET simulation at GRACEnet. Model AET estimates were more agreeable with ETIdaho at the LT manure location with a PBIAS of –14.3% and an NSE of 0.25; note that all research locations used identical weather files as they are near in proximity. Calibration endeavors only had a substantive impact on the LT manure location, NSE improved from 0.25 to 0.75, indicating good agreement. At GRACEnet, problems encountered with simulation of alfalfa growth require further attention to accurately reflect water use. In DNDC, the Penmen-Montieth equation is used to calculate daily potential evapotranspiration (PET); potential transpiration is then calculated based on crop demand for water and crop growth rate. If enough water is present and there is adequate supply of soil N, PET is met (i.e., AET = PET), otherwise water stress is incurred and crop growth is reduced. In ETIdaho, PET is calculated using the Penmen-Montieth equation and AET is calculated using the FAO-56 dual crop coefficient procedure that considers the effects of irrigation and precipitation surface wetting on evaporation [90]. While ETIdaho estimates are not direct field measurements, the AET estimation it employs has been shown to be accurate over a variety of crops in the western U.S. [33]. Guest et al. [89] reported that plans were in development to refine DNDC estimates of ET. In the present study, it was evident that problems encountered with DNDC simulated soil water and crop water use prevented their accurate simulation in a default capacity and were not acceptably ameliorated by calibration attempts. Additionally, though not considered in the present study, as the water balance underpredicted ET, leaching and/or runoff was likely overestimated.

3.5. DNDC Model Forecasting

The calibrated DNDC models were run for a subset of treatments at each location until 2050 to project management influence on long-term SOC stocks. Relative differences in SOC stock between “high” and “low” future emissions scenarios were 1.5%, 1%, and 1% at CoverCrop, GRACEnet, and LT manure respectively (Figure 6).

As an additional exercise, the simulation length was extended to year 2100 under the 52A treatment at LT manure which resulted in a more extensive relative difference (7.7%) between emission scenarios with the “high” scenario resulting in lower SOC (Figure S36).

The LT manure projection indicated no long-term advantage to applying manure annually or biennially when the net quantity of application was equal, supporting our empirical analysis. Soil organic C stocks of the 18A and 36B treatments best fit quadratic relationships (all R² > 0.85 and all p < 0.001), model critical points suggest an equilibrium will be reached at 63 Mg ha^–1 obtained by 2057. The 52A treatment also had a quadratic relationship (R² > 0.95, p < 0.001), the modeled critical point suggested an equilibrium of 130 Mg ha^–1 obtained by 2058. An additional simulation was carried out until year 2100 for the 52A treatment. This simulation indicated that an equilibrium level was not, in fact, realized as SOC contents were in excess of 200 Mg ha^–1 (Figure S36). This is incompatible with regional soils data which indicate that long-term (20–30 years records) manure application results in SOC levels between 80 and 120 Mg ha^–1 with only a handful of sites obtaining > 160 Mg ha^–1 SOC (Megan Satterwhite, IDA, personal communication). Considering all 2050 simulations, SOC was always within the levels of regional soil data (<160 g kg^–1); nevertheless, overestimation is possible as an equilibrium level was not explicitly obtained during the modeling period. In a comparable 45 years DNDC simulation of a ryegrass forage system, Khalil et al. [91] reported that SOC equilibrium was not obtained under several simulated management practices including additions of cattle and swine manure. Elsewhere in other long-term DNDC simulations, equilibrium states were suggested or obtained, but overall C input was much lower than in the present study [92,93]. Therefore, it is possible that when C inputs are very large the DNDC model projected SOC dynamics may erroneously favor sequestration which could inaccurately influence various NetZero C initiatives. The DNDC model may require integration of a “SOC ceiling” for variations in management, climate, and physical soil properties or a tuning of the process-based feedbacks that should regulate SOC levels. Additionally, when possible model projections should be supported by long term data records as is best practice.

At GRACEnet, fall manure applications resulted in 32% higher SOC stock relative to spring applications by 2050, likely resulting from unintentional imbalanced manure-C input which was exacerbated over time due to the cycled yearly model input (Table 2). Under synthetic fertilizer applications SOC was projected to decline over the next 30 years at a linear rate of 0.25 Mg ha^–1 year^–1at both LT manure and GRACEnet and 0.11 Mg ha^–1 year^–1 at CoverCrop except where triticale cover and minimum-till were utilized concurrently (Figure 6). Theoretically, synthetic fertilizer application should result in SOC accrual only if a net increase in C input is attained through increased crop growth. At CoverCrop, the addition of triticale cover would provide additional C input, yet DNDC projects this is not sufficient to positively effect SOC levels under a more disruptive tillage regime. Here, the MT-NM-T treatment accumulated SOC following a quadratic model (y = –34302 + 34x –0.0082x², R² = 0.9, p < 0.001), while SOC of the CT-NM-T treatment declined linearly at 0.08 Mg ha^–1 year^–1. Under manure applied treatments at CoverCrop, projected SOC similarly followed strong quadratic relationships (all R² > 0.95 and all p < 0.001); evaluation of model critical points suggested new SOC equilibriums would be met between 2055 and 2056 at SOC levels between 94 and 121 Mg ha^–1. By 2050, SOC stock of conventionally tilled treatments receiving manure were 17% higher on average than their minimum-till counterparts. Contrasting this, when no manure was applied minimum-till treatments were 8% higher on average than conventionally tilled treatments. Continuing comparison of CoverCrop, projected SOC of manure applied treatments employing triticale cover were 8% higher on average than corresponding fallow treatments by 2050; when no manure was applied, treatments utilizing triticale cover were 15% higher on average. Thus, the DNDC projections indicate a long-term benefit of triticale cover regardless of tillage or manure application, and a variable effect of tillage practice dependent on manure application. Our projections on tillage support the meta-analysis by Gross and Glaser ([94], Under Review) and the reasoning of Baker et al. [95] that suggest SOC accumulation appears to favor conventional tillage management, at least to incorporate manure, when it is applied. The effect of adding triticale cover can be calculated as the absolute difference in the rate of SOC change in no-manure treatments as projected by DNDC; these estimates (conventional tilled = 0.05 Mg ha^–1 year^–1 and minimum-till = 0.27 Mg ha^–1 year^–1) are comparable to the annual rate of change (0.32 ± 0.08 Mg ha^–1 year^–1) reported in a recent meta-analysis of cover cropping [66]. Poeplau and Don [66] included only those studies where winter cover was not harvested which may explain our smaller estimates.

4. Conclusions

Empirical measurements indicated a commercial dairy forage rotation of corn-barley-alfalfa^x3 accumulated SOC (0.6 ± 0.5 Mg ha^–1 year^–1) while a commercial rotation of wheat–potato–barley–sugarbeet had not significantly changed in 8 years excluding the effects of manure. Manure application timing (fall vs. spring) and frequency (1X annual vs. 2X biennial) did not significantly affect SOC accrual (p = 0.41 and p = 0.75) when the total amount applied was equivalent. These findings better inform crop producers in semi-arid climates utilizing manure resources and may also influence attempts to sequester C in similar regions. The DNDC model simulated observed values of SOC and biomass C acceptably in a “default” capacity throughout the duration of observed data, but problems were encountered with simulation of soil water contents and AET. Difficulties in simulation of soil water suggests further evaluation of modelled GHG emissions (CH₄, CO₂, N₂O) in semi-arid regions. Forecasting of calibrated models to 2050 suggests triticale cover accrues SOC (0.05–0.27 Mg ha^–1 year^–1), and that conventional (disk/chisel) tillage is favored over minimum-till (relative SOC difference = 17%) when manure is applied. Forecasting SOC of manure applied treatments suggested quadratic relationships (all R² > 0.85 and all p < 0.0001), however when extended to year 2100 no equilibrium was realized. Modeling results provide a reference for users of the “default” DNDC model in similar regions and indicate that overestimation of SOC sequestration potential is possible.