Modeling Hairy Vetch and Cereal Rye Cover Crop Decomposition and Nitrogen Release

: Empirical models could help us to understand the process of plant residue decomposition and nutrient release into the soil. The objective of this study was to determine an appropriate model to describe the decomposition of hairy vetch ( Vicia villosa Roth) and cereal rye ( Secale cereale L.) cover crop (CC) residue and nitrogen (N) release. Data pertaining to above and belowground CC residue mass loss and N release for up to 2633 cumulative decomposition degree days (112 d) after litterbag installation were obtained from two cropping system experiments, a 1-yr study conducted in 2015 and a 2-yr study during 2017 to 2018 in the humid subtropical environment of southern IL, USA. Six exponential and two hyperbolic models were fit to percent mass and N remaining data to find the one with minimum Akaike information criterion (AIC) and residual sum of squares. Modified three-parameter single exponential and two- or three-parameter hyperbolic models best met the assumed criteria of selection for above and belowground CC residue, respectively. Fitting a double exponential model to combined data for percent mass and N remaining identified two mass and N pools, a fast and a slow pool with different rate constants. A five-parameter double exponential with an asymptote met the preset criteria and passed all tests for normally distributed population, constant variance, and independence of residuals at α = 0.05 when fit to combined data of hairy vetch shoot mass and N remaining. However, a two-parameter hyperbolic and three-parameter asymptotic hyperbolic model provided the best fit to a combined data of cereal rye shoot mass and N remaining, respectively. Both hyperbolic decay models showed a good fit for belowground mass decomposition and N release for both CCs. Cereal rye had a poorer fit than hairy vetch for mass and N remaining of both above and belowground mass. The best-selected decay models can be used to estimate the decomposition and N release rates of hairy vetch and cereal rye above and belowground residue in a similar environment. writing—review


Introduction
Cover crop (CC) residue is a source of soil organic matter, and its degradation is critical to subsequent crop productivity. Residue decomposition determines the soil nutrient pool and regulates nutrient release in soil [1] through depolymerization of fibers and hydrolysis of sugars, mostly via heterotrophic soil microorganisms [2]. External environmental conditions as modified by climate, agronomic practices such as tillage, cropping rotation, fertilizer, and irrigation, and plant residue quality may impact the rate of decomposition and nutrient release. Inherent properties of the residue such as carbon-nitrogen (CN) ratio, fiber fractions, and lignin concentration can greatly affect the litter decomposition and nutrient cycling [3][4][5]. These properties differ between C3-and C4-derived soil organic matter [6] and between grass and legume residue [7], which may impact decomposition and nutrient release kinetics, indicating the possibility of the usefulness of the different approaches for modeling those kinetics. The choice of approach also depends on the desired degree of analytical simplicity, predictive power, and generality [8]. Knowledge of decay mechanism and use of a suitable model specific to the substrate quality can provide valuable information for CC management, which is mostly lacking in comparative studies where a single model is chosen for a variety of crops to determine decay rate constants and half-lives. There is a lack of uniformity in using decay models for decomposition and mineralization studies for similar substrates, which vary from simple one parametric single exponential first-order models to the complex multiparametric consecutive exponential models.
A first-order single exponential decay model (y = ae −bx , where y is the mass of substrate at time x, b is the rate constant, and e is the base of the natural logarithms (2.71828)) has been widely used for nutrient mineralization, residue decomposition, and plant population studies [5,7,[9][10][11][12]. It was applied for modeling litter decomposition for numerous grasses and legumes [7] and fine litter decomposition of forest soil [12]. Ruffo and Bollero [13] and Sievers and Cook [5] used this model with an asymptote (y = ae −bx + yₒ, where yₒ is an asymptote) for cereal rye (Secale cereale L.) and hairy vetch (Vicia villosa Roth) decomposition and nutrient release. Polglase et al. [14] used a single exponential model for P mineralization from soil organic matter in a pine forest, whereas Fernández et al. [15] used it for modeling C mineralization in soils after wildfires in Spain. The strength of this model is that it produces a single rate constant, which can be used directly to compare decay rates from different treatments. However, it does not accurately describe decomposition or mineralization kinetics, where rate constants vary with time due to rapid loss or an extended lag phase in early decomposition [9,12,16]. The CC-derived labile fraction of soil organic matter composed of light (low specific density or mineral-free) and heavy (high specific density or mineral-bound) fractions tends to follow a different kinetic model in describing the decomposition and mineralization [17,18].
The first-order double exponential model with two rate constants (y = ae −bx + ce −dx , where b and d are the rate constants) separates organic matter into a soluble fraction (e.g., sucrose) or fast pool and cell-wall (e.g., detergent fibers) or slow pool [19] fraction. It was reported to have improved goodness of fit of single exponential models for residue decomposition and nutrient release mechanisms [9,18,20]. Berndt [9] suggested this model over the single exponential model when comparing kinetic parameters of decay of C remaining for hybrid bermudagrass (Cynodon dactylon (L.) Pers. × Cynodon transvaalensis Burtt-Davy) thatch. Wang et al. [20] predicted temperature-and moisture-dependent rate constants for soil N mineralization with a modified double exponential model under standard temperature (35 °C) and moisture conditions (55% water holding capacity). Dhakal et al. [11] reported that the double exponential model described alfalfa (Medicago sativa L.) population decline in a semiarid environment, which resulted in the highest adjusted R 2 (0.94 to 0.97) and the lowest standard error of estimate (SEE) for the upright-type alfalfa cultivars. Fernandez et al. [15] and Camargo et al. [21] reported this model as fitting mineralization data better than a single exponential model. Although all models generated an R 2 greater than 0.98 for in vitro mineralization of C, the double exponential model could not fit some of the samples, whereas exponential with a linear combination (y = ae −bx + cx + yₒ, where c is the slope of the linear function) yielded superior results to the double exponential, exponential plus an asymptote, and hyperbolic model [y = ab/(b + x)] [17]. Dendooven et al. [22] reported poor fit of the double exponential function in fitting N mineralization data to characterize active and recalcitrant organic N pools derived from sugar-beet (Beta vulgaris L.) and bean (Phaseolus vulgaris L.) residue.
Besides exponential models, a hyperbolic function was recommended to minimize standard errors, compared to the first-order exponential model in fitting the N mineralization data [16]. Decay and N release of cereal rye and hairy vetch residue has been well described using the hyperbolic model when compared to linear and first-order models [23]. In contrast, Berndt [9] reported poor fit statistics for the two-parametric hyperbolic decay function, relative to exponential models. This indicates the need to test various empirical models, specific to the plant species. Since mass loss and N release from cereal rye and hairy vetch residue have been studied using a variety of empirical models [5,13,23], performances of those functions have not been compared yet to suggest the best fit model.
The current study provides an overview of performances of the commonly used empirical models in CC decomposition and N mineralization studies. The objective of this research was to examine eight mathematical models to test their statistical significance in explaining cereal rye and hairy vetch decay and N mineralization in a sub-humid environment. Mass and N remaining of CC residue were fitted with six exponential and two hyperbolic models, and statistical parameters were compared for those models. An empirical model with the highest adjusted R 2 and lowest residual sum of squares and Akaike information criterion (AIC; [24]) may be considered best for decomposition studies.

Experiment Locations
Data from two experiments (Experiment 1, [5]; Experiment 2, [25]), comprised of two different intervals, were used to carry out this study. Both studies were conducted at the Agronomy Research Center (ARC, 37.7029 N, −89.2403 W and 38.185 N, −89.4592 W, respectively) at Southern Illinois University, Carbondale, IL. Soil series at both locations were Hosmer silt-loam (fine, silty, mixed, active, mesic Oxyaquic Fraguidalfs). Research design, treatments, site soil properties, and weather conditions were described in greater detail by Sievers and Cook [5] and Singh et al. [25]. Exp. 1 location received nearly 540 mm cumulative rainfall (> 80% within 67 d after beginning the trial) during the study period in 2015, while 482-and 414-mm cumulative rainfall were recorded in 2017 and 2018 in Exp. 2, respectively. The Exp. 2 location received > 50% of the total rainfall within 12 d in 2017, whereas ~ 70% of the cumulative rainfall was received within the first two months in 2018. In Exp. 1, the maximum average daily temperature recorded was 34.9 °C on 6 July 2015 and the minimum was 5.

Sampling and Data Collection
In Exp. 1, CC biomass was collected in spring 2015 from two locations: cereal rye from Kuehn research farm (37.933548 N, −89.244436 W) of the Southern Illinois University, Carbondale, IL; whereas hairy vetch was obtained from the ARC. The experiment was overlaid as a portion of a larger trial where only the no-till subsystem plots were used for the decomposition study. Hairy vetch and cereal rye were terminated on 15 and 23 April 2015, respectively, using glyphosate [N-(phosphonomethyl) glycine] at 1.26 kg ae ha −1 . Aboveground biomass was clipped 3 to 4 d after burndown and air-dried for 1 wk to be used as a residue in the litterbags, whereas samples clipped before herbicide application were oven-dried to determine in situ CC production and C to N ratio. Belowground biomass was collected by taking intact root cores of 5 cm in diameter and 15 cm in length. Ten grams of air-dried aboveground CC biomass was placed into traditionally made 20 × 20 cm nylon mesh bags, with 5 mm mesh on the upper side and 2 mm mesh on the bottom. The aboveground biomass used was equivalent to 1277 and 2203 kg ha −1 dry mass and 14.7 and 92.3 kg N ha −1 with a C to N ratio of 34.7 and 9.5 for cereal rye and hairy vetch, respectively, which was much greater than the actual in situ productivity, i.e., 578 and 791 kg ha −1 dry mass and 13.5 and 21.4 kg N ha −1 for cereal rye and hairy vetch, respectively [5]. A total of 14 litterbags were installed in each notill sub-plot under soybean [Glycine max (L.) Merr.] and corn (Zea mays L.) main plot, which were rotated every year in four replicates, giving a total number of 112 litterbags (56 cereal rye + 56 hairy vetch). Litterbags were placed on the ground on 5 May 2015, and biomass samples were collected on the same day for 'week 0′ sampling. After that, two litterbags per plot were collected at 2, 4, 6, 8, 12, and 16 wk after litterbag installation. Corn and soybean were planted on the 4th and 12th of June 2015, respectively, and the growth stage of the crops was noted at the time of litterbag collection. The intact root cores remained in plastic liners, wrapped with 16 mesh size at open ends were used for belowground study. The same number of cores as the litterbags were inserted 5 cm below the soil surface and collected on the same day as mentioned for the litterbags. Collected soil cores were washed to retrieve belowground biomass then oven-dried and ground to analyze for its constituents. Corn was fertilized with N at the late vegetative stage [5]. Aboveground biomass and collected residue samples for each CC were dried, ground to pass through a 1-mm screen, and analyzed for lignin, acid detergent fiber, neutral detergent fiber, and C and N concentration. The complete field events were described by Sievers and Cook [5] in greater detail.
Similarly, Exp. 2 was laid out in a completely randomized design with three replicates, overlapped on an ongoing tillage study established in fall 2013. The experiment consisted of two CCs (cereal rye and hairy vetch) under two tillage systems (reduced tillage and no-tillage), giving a total of 12 plots. Corn and soybean were grain crops in 2017 and 2018, respectively. Cover crops were killed by spraying a mixture of glyphosate at 1.26 kg ae ha −1 , 2,4-D (Dichlorophenoxyacetic acid) at 0.80 kg ae ha -1 , and ammonium sulfate (2.5% v/v) on 13 and 28 April in 2017 and 2018, respectively. Aboveground CC biomass was collected 2 to 3 d after spraying, washed, and air-dried for 1 to 3 d. Fifty grams of air-dried samples were packed into the litterbags of the mesh size described above for Exp. 1. The packed CC biomass was equivalent to 2206 to 2743 kg ha −1 for cereal rye and 1198 to 1624 kg ha −1 for hairy vetch dry biomass. A total of 132 litterbags were used in each year. At "wk 0" all litterbags were placed on the ground to simulate installation and returned to the lab for further analysis. Six days after "wk 0" in 2017 and 8 d after "wk 0" in 2018, tillage operation was performed at reduced treatment plots to plant the main season crop. After that, litterbags in tillage treatment were placed in a vertical opening down to 15 cm into the soil. The no-till treatment had litterbags on the soil surface throughout the study in both years. In 2017, litterbags were installed on 19 April (week 0) and then one litterbag per plot was collected weekly for 10 wks, whereas in 2018, they were placed on 2 May and collected at 1, 2, 3, 4, 5, 6, 8, 10, 12, and 14 wks after installation. Sample preparation, grinding, and lab analysis was similar to the method described above.
For both experiments, the percentage of ash-free mass remaining (MR, %) and N remaining (NR, %) at a given time was calculated using the formula: MR or NR = (Xt/Xo) × 100 (1) where X was the mass or N at a given time t (decomposition degree days, DDD), and Xo was the initial CC mass or N mass at week 0. To standardize time after litterbag installation based on daily air temperature and DDD, it was calculated as follows [26]: where TMax and TMin are the daily maximum and minimum air temperature, respectively, TBase is the base temperature for the CC decomposition, which was considered 0°C [26]. When TMax or TMin were less than TBase, the TMax and TMin computed were equal to TBase. For the days when TMax was greater than 30°C, the TMax was changed to 30°C.

Comparison of Empirical Models
Eight non-linear models were fitted to the percent mass and N remaining vs. accumulated DDD. One of the first-order decay models tested was a two-parameter single exponential decay model by Stanford and Smith [27], which captures the gradually slowing absolute rate of mass loss over time at a constant temperature and moisture [28]: where a is the y-intercept or numeric constant to satisfy the model, b is the relative decay rate or proportionality constant, and x is an independent variable or time. Howard and Howard [29] and Wieder and Lang [30] added an asymptote (yₒ) to capture the resistant litter fraction (Equation (4)).
A modified three-parameter single exponential decay model [9,11] has also been used to compare with other exponential models as provided by Systat Software [31] (Equation (5)).
where c is the numeric constant. Single exponential models have been criticized for not representing the transition from rapid to slow decomposition, whereas the double exponential model with two single exponential components, which consists of two decay or mineralization rate constants, is reported to be a better alternative [9,19,32]. A four-parameter double exponential model can be written as Equation (6) [33].
where a and c are the constants and b and d are the rates of decay of available (light) and resistant (heavy) fractions of residue, respectively. An asymptote can be added to Equation (6) to further catchup the resistant fraction of the residue. The five-parameter function used was [31] y = ae −bx + ce −dx + yₒ A double-pool model was reported in yielding significantly smaller root mean square errors in which one pool was assumed to mineralize exponentially and the second pool worked according to zero-order kinetics [10,34] for modeling the flush of N mineralization caused by drying and rewetting soils.
where c is also the rate constant for the mineralization of the slow pool fraction of the residue. Besides exponential decay models, hyperbolic equations were also found effective in explaining N mineralization in soils [16]. The two-parameter hyperbolic decay model tested in our study was where b is the rate constant. The three-parameter hyperbolic model with asymptote was also used for comparison [31].
The data were subjected to a Lavene's test and Shapiro Wilk test for variance and normality of data at α = 0.05, respectively, using PROC NLIN in SAS 9.4 (SAS Institute, Cary, NC). In addition, partial residual plots for mass and N remaining against time were visually analyzed to confirm the non-linear pattern of data. Then, models were fitted for percent mass remaining and N remaining for each of the studies and CCs using SigmaPlot 14.0 [31]. For Exp. 1, models were fitted for aboveground and belowground root biomass. Data from two tillage treatments were combined within each CC for Exp. 2 for both study years. The iterative method adopted in SigmaPlot was based on the Marquart-Levenberg algorithm [35] for all non-linear models. Models were compared based on normality, the constant variance test [31,36], and the Durbin-Watson test [37] to detect positive or negative autocorrelation of residuals. These tests were conducted at α = 0.05, where the models were assumed to be passed or failed based on a given standard criterion. Test statistics such as adjusted R 2 , standard error of estimate (SEE), residual mean squares (RMS), predicted residual error sum of squares (PRESS), and Akaike information criterion was also used for model comparison. Model fitting excluded influential outliers using Leverage and Cook's D.
The Akaike information criterion is good for model selection; however, with the increase in complexity of the model, such as from single to multiple exponential functions, the AIC may fail to select the best model because the criterion assumes that the true model is among the candidate pool, in a condition that none of the models are representing a complete set of data. To solve the problem, the PRESS statistic has often been used for cross-validation of models, which uses a predicted set of samples to provide an unbiased evaluation of predictability of the model [38]. Models passed normality, variance test, and residual test with the highest adjusted R 2 , and the lowest SEE, RMS, PRESS, and AIC values were considered the best fit for hairy vetch and cereal rye CC decomposition and nutrient mineralization. Model parameters were estimated for each species, year, and study. Regression plots were obtained from SigmaPlot 14.0 [31].

Modeling Percent Mass Remaining
Statistical values and parameters of eight different non-linear models explaining the percent mass remaining of CC residues are given in Table 1-4. All models were valid in predicting mass loss (P < 0.001). For hairy vetch aboveground residue in Exp. 1, all models had an R 2 value of 0.97 except for the two-parameter single exponential decay model (0.91) ( Table 1). The modified three-parameter single exponential function had the lowest RMS, SEE, PRESS, and AIC values. Although fiveparameter double exponential and hyperbolic decay models had SEE comparable to the modified single exponential model, these models failed the tests for normality and independence of residuals (Durbin-Watson), whereas the modified single exponential model passed those tests, including the constant variance of errors. The four-parameter double exponential model produced greater R 2 and lower SEE and PRESS statistics, while the four-parameter single exponential with linear combination resulted in a comparable R 2 and SEE to the double exponential model, and lower RMS and AIC for cereal rye aboveground biomass (Table 1). However, the latter failed the normality test and the test for independence of residuals. The double exponential model passed all those test criteria and appeared to be a promising model for aboveground cereal rye residue decomposition. The fiveparameter double exponential model with an asymptote had non-significant rate constants, especially for the resistant fraction of the cereal rye residue.
The model that reduced AIC and PRESS statistics in Table 2 was the two-parameter hyperbolic decay for hairy vetch belowground biomass decay. Nevertheless, the four-parameter double exponential model best minimized the RMS and SEE. However, double exponential models had a non-significant rate constant (P > 0.01) for slow pool fraction and failed assumption of normally distributed population. Belowground mass remaining for cereal rye was also explained better by the two-parameter hyperbolic decay model, which best minimized RMS and SEE with the highest adjusted R 2 and the lowest AIC and PRESS statistic ( Table 2). The model also satisfied the assumption of normally distributed population, constant variance, and independence of residuals. Double exponential models had at least one of the parameters as non-significant in predicting the hairy vetch and cereal rye mass decomposition. Table 3 compares the models for aboveground biomass decomposition for hairy vetch and cereal rye from Exp. 2. Model performances did not vary between tillage treatment when the parameters were compared (P > 0.05). Across the tillage treatments, the five-parameter double exponential model with an asymptote best minimized the RMS and SEE and gave the lowest AIC and PRESS, relative to other close models for hairy vetch decomposition. It also passed the tests for normality, variance, and independence of residuals. All models generated the same adjusted R 2 value (0.90) except for a simple single exponential function (0.86) for hairy vetch ( Table 3). None of the models passed all tests for normality, constant variance, and independence of residuals for cereal rye percent mass remaining. The five-parameter double exponential model passed the tests for normality and constant variance at α = 0.05 but produced non-significant estimates of parameters. The single exponential with an asymptote yielded the lowest PRESS and AIC and passed tests for normality and independence of residuals at α = 0.05 with significant model parameters. This single exponential model was found to be best for cereal rye in Exp. 2.  2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where pass or fail assumptions were made at α ≤ 0.05. 9 Not-a-number notation for non-finite residuals. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level. ; 2, y = ae −bx + yₒ; 3, y = ae b/(c + x) ; 4, y = ae −bx + ce −dx ; 5, y = ae −bx + ce −dx + yₒ; 6, y = ae −bx + cx + yₒ; 7, y = ab/(b + x); 8, y = ab/(b + x) + yₒ. 2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where pass or fail assumptions were made at α ≤ 0.05. 9 Not-a-number notation for non-finite residuals. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level. ; 2, y = ae −bx + yₒ; 3, y = ae b/(c + x) ; 4, y = ae −bx + ce −dx ; 5, y = ae −bx + ce −dx + yₒ; 6, y = ae −bx + cx + yₒ; 7, y = ab/(b + x); 8, y = ab/(b + x) + yₒ. 2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where pass or fail assumptions were made at α ≤ 0.05. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level.  2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where pass or fail assumptions were made at α ≤ 0.05.
We compared the fitness of exponential and hyperbolic functions to combined data from two studies ( Table 4) as portrayed by Figure 1. The shape of the decay models followed a pattern of rapid mass loss from day 0 to nearly 1000 accumulated DDD and a slow rate of decomposition afterward ( Figure 1). All models produced very high adjusted R 2 and low SEE values except a two-parameter single exponential model for both CC residues. None of the models passed all three statistical tests viz. test for normality, constant variance, and independence of residuals for both hairy vetch and cereal rye. Results showed better fit with five-parameter double exponential plus an asymptote than the single exponential and hyperbolic models for hairy vetch CC decomposition. However, the twoparameter hyperbolic model also produced standard errors and AIC values close to the fiveparameter double exponential model in minimizing RMS, SEE, and AIC. Despite that, the choice between the five-parameter double exponential and two-parameter hyperbolic model would suggest the former model as the best fit with significant heteroskedasticity (Table 4). In contrast to the exponential models, the two-parameter hyperbolic model seemed to have the best fit for the cereal rye percent mass remaining data, as the SEE, RMS, PRESS, and AIC appeared lower than or equal to exponential and three-parameter hyperbolic models. This model also passed the test for the constant variance of the errors. Overall, the five-parameter double exponential model with an asymptote appeared suitable for hairy vetch decomposition modeling, whereas cereal rye had inconsistent results for individual small datasets and the two-parameter hyperbolic model for the combined data. Exponential and hyperbolic decay models explaining percent mass remaining of hairy vetch and cereal rye cover crop aboveground residue against cumulative decomposition degree days at 0 °C base temperature. (i) Two-parameter single exponential, (ii) three-parameter single exponential, (iii) modified three-parameter single exponential, (iv) four-parameter double exponential, (v) fiveparameter double exponential, (vi) four-parameter single exponent with linear combination, (vii) two-parameter hyperbolic, and (viii) three-parameter hyperbolic. The upper equation represents hairy vetch and the lower cereal rye. Data were pooled from Exp. 1 and 2. All models were significant at P < 0.0001.

Modeling Percent Nitrogen Remaining
Nitrogen released from above-and belowground CC residue was non-linear with accumulated DDD (P < 0.001). The results followed a similar pattern to those of the percent mass remaining of the residue. Table 5-8 describes the parameter estimates and test statistics for the six exponential and two hyperbolic decay models. Table 5 shows results from Exp. 1 for percent N remaining of aboveground CC residue. The adjusted R 2 of the models was near perfect (>0.96) while for cereal rye it ranged from 0.67 to 0.70 when fit to percent N remaining data ( Table 5). All models passed the constant variance of residuals test for both CCs except for the three-parameter single exponential model with an asymptote. The modified three-parameter single exponential model appeared to have the best fit for the hairy vetch percent N remaining, which minimized RMS and SEE and lowered the PRESS and AIC statistics, relative to other decay models. For the aboveground cereal rye percent N remaining, the four-parameter single exponential model with the linear combination had the highest adjusted R 2 (0.70) and the lowest RMS, SEE, PRESS, and AIC values in relation to other exponential and hyperbolic models (Table 5), but the rate constant was not significant. This means the model cannot explain the N release rates. Thus, the three-parameter single exponential model was chosen based on relatively smaller SEE, RMS, and AIC and greater adjusted R 2 . This model also passed assumptions for the normal population, constant variance, and independence of residuals. ; 2, y = ae −bx + yₒ; 3, y = ae b/(c + x) ; 4, y = ae −bx + ce −dx ; 5, y = ae −bx + ce −dx + yₒ; 6, y = ae −bx + cx + yₒ; 7, y = ab/(b + x); 8, y = ab/(b + x) + yₒ. 2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model, 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where pass or fail assumptions were made at α ≤ 0.05. 9 Not-a-number notation for non-finite residuals. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level.  15.56 1 Decay models: 1, y = ae −bx ; 2, y = ae −bx + yₒ; 3, y = ae b/(c + x) ; 4, y = ae −bx + ce −dx ; 5, y = ae −bx + ce −dx + yₒ; 6, y = ae −bx + cx + yₒ; 7, y = ab/(b + x); 8, y = ab/(b + x) + yₒ. 2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where a pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where a pass or fail assumptions were made at α ≤ 0.05. 9 Not-a-number notation for non-finite residuals. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level.  1 Decay models: 1, y = ae −bx ; 2, y = ae −bx + yₒ; 3, y = ae b/(c + x) ; 4, y = ae −bx + ce −dx ; 5, y = ae −bx + ce −dx + yₒ; 6, y = ae −bx + cx + yₒ; 7, y = ab/(b + x); 8, y = ab/(b + x) + yₒ 2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where a pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where a pass or fail assumptions were made at α ≤ 0.05. * t-test significant at α ≤ 0.01, ** at α ≤ 0.001, and ns, not significant at the α = 0.01 level.  2 Residual mean square of the model. 3 Standard error of estimate of the model parameters. 4 Predicted residual sum of squares estimate of the model. 5 Akaike information criterion value of the model. 6 Shapiro-Wilk test for normality of the data where pass and fail assumptions were made at α ≤ 0.05. 7 Constant variance test using Spearman rank correlation where a pass or fail assumptions were made at α ≤ 0.05. 8 Durbin-Watson test of independence of residuals where a pass or fail assumptions were made at α ≤ 0.05. Table 6 shows fit parameters and statistics for below-ground CC residue from Exp. 1. Four-and five-parameter exponential models had non-significant decay rate constants for the slow pool of the residue and had relatively higher RMS and AIC than hyperbolic functions. The two-parameter hyperbolic model produced high adjusted R 2 and minimized RMS, SEE, PRESS, and AIC for hairyvetch N remaining data when compared to exponential and three-parameter hyperbolic decay models. For cereal rye N remaining, three-parameter hyperbolic decay function with an asymptote fitted best in minimizing RMS, SEE, and AIC, while the model also passed the assumption of normality, constant variance, and independence of residuals.
All models fitted to percent N remaining data from Exp. 2 failed the normality test (Table 7); however, some of the models passed the constant variance and independence of residual tests. The five-parameter double exponential model with an asymptote fitted to the percent N remaining data for hairy vetch best minimized the RMS, SEE, PRESS, and AIC with the greatest adjusted R 2 value and significant rate constants. The model also passed a test for constant variance and independence of residuals. The model that best minimized the RMS and SEE for the percent N remaining of cereal rye was the modified three-parameter single exponential. The model produced the greatest R 2 and had the lowest AIC value. However, the tests for normality, variance, and residuals were not satisfied by any of the models for cereal rye.
Similar to the individual studies, the five-parameter double exponential model fitted best for the combined dataset with very high adjusted R 2 (0.94) and the lowest SEE, RMS, PRESS, and AIC values for hairy vetch N remaining (Table 8). All models failed the test for normality for both CCs. The fiveparameter double exponential function passed a test for constant variance and independence of residuals. For cereal rye percent N remaining, the three-parameter hyperbolic model with an asymptotic best minimized the RMS and SEE and had the lowest PRESS and AIC values ( Table 8). None of the models could satisfy the assumption of normality, constant variance, and independence of residuals. The double exponential model also produced high adjusted R 2 and minimized the RMS and SEE, but had non-significant rate constants for percent N remaining of cereal rye residue. The modified three-parameter single exponential model was equally as good as the hyperbolic decay model. Figure 2 visualizes the pattern of those selected exponential decay models where more than 80% of the hairy vetch N was mineralized within the first 600 accumulated DDD and nearly 50% cereal rye N was mineralized within 1000 accumulated DDD from the period of total 2700 DDD. It indicates that there were two phases of N release into the soil: progressive and lag. The rapid N release rate during the progressive phase took less than 25% of the total DDD for hairy vetch and less than 40% of the total DDD for cereal rye CC residue. Overall, hairy vetch N dynamics clearly followed double exponential function, while the cereal rye exhibited mostly the hyperbolic decay function for mass and N remaining. Cereal rye produced higher residuals than hairy vetch because of more spread of the data, especially during the initial decomposition period. Exponential and hyperbolic decay models explaining percent nitrogen remaining of hairy vetch and cereal rye cover crop aboveground residue against cumulative decomposition degree days at 0 °C base temperature. (i) Two-parameter single exponential, (ii) three-parameter single exponential, (iii) modified three-parameter single exponential, (iv) four-parameter double exponential, (v) five-parameter double exponential, (vi) four-parameter single exponent with linear combination, (vii) two-parameter hyperbolic, and (viii) three-parameter hyperbolic. The upper equation represents hairy vetch and the lower cereal rye. Data were pooled from Exp. 1 and 2. All models were significant at P < 0.0001.

Discussion
Results showed a stronger relationship between percent mass or N remaining (high adjusted R 2 , Table 1-8) of hairy vetch CC residue and cumulative DDD than that of cereal rye residue (Figures 1  and 2). The decomposition rate constants and asymptotes were dissimilar for these two CCs (Table  1-8). Although Exp. 2 had two tillage systems (reduced vs. no-tillage), there was no effect on decomposition rate constant and N release in 2017 and 2018 [25], indicating that the management factor had lesser impact than the substrate quality factor. The inherent plant chemical constituents determine the rate of decomposition and N mineralization into the soil [4]. In both CC decomposition experiments, the CN ratio was greater for cereal rye (35:1 and 24:1 for Exp. 1 and 2, respectively) than hairy vetch (10:1, 9:1 for Exp. 1 and 2, respectively) [5,25]. Greater N concentration and less fiber content in hairy vetch may have accelerated the decomposition and N release in the early days or with less accumulated DDD, a relationship reported by Ruffo and Bollero [13] and Otte et al. [39]. Nitrogen and soluble carbohydrate fraction of the plant residue enhanced the microbial growth efficiency and improved the bermudagrass residue decomposition in FL, USA [9]. In both experiments, the residuals were higher for the cereal rye mass remaining and N release data than for hairy vetch. The lack of fit was due to large variability in initial mass and N content of cereal rye and immobilization of the N in the first 4 wks [5]. Rufo and Bollero [13] also reported the lack of fit for cereal rye when compared to hairy vetch. The presence of high neutral detergent fiber and lignin concentration [5] likely influenced the tensile strength of leaves and decomposability of cereal rye, which provide mechanical and chemical defense against microbial and chemical degradation [40]. Cornelissen et al. [41] reported that leaf tensile strength was related to litter decomposition for C3 grasses and also suggested that the complex leaf base content of the grass species can result in high variation in mass and N loss.
The greater adjusted R 2 values and lower standard errors with the double exponential model with two rate constants when compared to a single exponential with or without an asymptote, especially for hairy vetch CC residue, indicated two residue pools, i.e., fast and slow decomposing fractions. This could be due to the presence of labile versus and recalcitrant fractions of the plant materials. However, in some cases, the addition of asymptotic or linear or another exponential component to the simple single exponential model resulted in non-finite residuals and failed to crossvalidate the model, mainly for cereal rye for individual datasets from Exp. 1 and 2. In such cases, the slopes (or rate constants) for the double exponential model or exponential model with the linear combination had poor predictability at α ≤ 0.01. This indicates that the sample size was too small to have a sufficient sampling frame for the model, which gave undue weight to the initial data points for the short study period. We noticed that the issue of non-finite residuals was eliminated with combined data from Exp. 1 and 2, which extended the x-axis from 14 wk (cumulative DDD = 2382) to 16 wk (cumulative DDD = 2633) and increased the number of XY pairs. That was the reason Otte et al. [39] suggested using the simple asymptotic model rather than a double exponential model for cereal rye biomass decomposition. However, an asymptotic model may not always give the best results when compared to non-asymptotic models. The plant chemical constituents and environmental conditions may also affect the model performances [40,42]. Both of our experiments received more than half of the total cumulative rain in the first few weeks and the volumetric soil water content was near the field capacity, which might help accelerate the residue degradation as most of the soil microbes are highly active and thrive under moist warm conditions [42]. Hairy vetch aboveground biomass was well represented by a double exponential model with an asymptote with high R 2 and minimal residual errors and AIC with shorter fast-pool turnover time (~25% of the total accumulated DDD), compared to that of cereal rye (~40% of the total accumulated DDD). Juma et al. [16] suggested that a model of N cycling in soil should have at least three to four compartments to account for different sources of N.
Only hyperbolic models produced relatively low RMS, SEE, AIC, and PRESS statistics for the belowground hairy vetch and cereal rye percent mass and N remaining, while exponential models showed poorer fits. The fast-pool turnover time of belowground residue mass was shorter than that of aboveground biomass owing to its immediate proximity to soil biotic and abiotic environments [43][44][45]. Similar to aboveground biomass, there was a difference in mass loss and N release from belowground mass between hairy vetch and cereal rye. However, the pattern of decaying was different for belowground biomass even though the CN ratio of root mass was in line with aboveground litter, which was 17:1 for hairy vetch and 40:1 for cereal rye [5]. Sun et al. [43] found distinct traits other than the CN ratio that controls root decomposition where N was not the major driver, but C compounds associated with mycorrhiza were the major factors. They reported a poorer fit of the double exponential model for the root mass remaining data. The parameter estimates vary more with the model fitting procedure for exponential equations than the hyperbolic models [16]. The reason is that the first-order equation assumes that the rate of mineralization is proportional to the initial size (MR0 or NR0) of the mineralizable pool [27], which might be affected by environmental and management factors and generate more variance. In our study, the belowground initial mineralizable pool was considerably smaller than the aboveground pool, and a significant portion of root mass remained undecomposed after 550 DDD [5]; thus, exponential models produced an ambiguous estimate of parameters. The evidence of a wide range of parameter estimates of the firstorder equation was also reported by Nicolardot et al. [46] and Talpaz et al. [47].
The asymptotic approach used by Sievers and Cook [5] and Singh et al. [25] in Exp. 1 and 2, respectively, to estimate decomposition and N mineralization rate had higher RMS and SEE than the double exponential and hyperbolic models for both crops. The point of contention was that the single exponential asymptotic model mostly departed from the assumption of normally distributed population, constant variance, and independence of residuals. However, the lignin content of the plant materials leaves a certain amount of highly recalcitrant compound in the late stages of decomposition, despite having a high rate of decomposition, which can typically be addressed by the asymptotic form of equations [43]. However, use of an asymptote in a single exponential model for a short study period may limit the first-order decomposable pool to a smaller fraction and might give undue weight to the slowly decomposing fraction than the relatively long period. It is difficult to infer the persistence of the recalcitrant fraction and conversion of residue into the soil organic matter based on mass loss and N remaining data [48,49]. The addition of C dynamics and microbial growth rate factors in the study would be a step towards the empirical calculation of soil organic matter conversion. However, the current study was able to detect the statistical precision and fitness of the commonly used models in decomposition and N dynamics studies. Finding control mechanisms that fitted hyperbolic function to the root decomposition data would extend the prospect of this study.

Conclusions
The double exponential model with an asymptote can be used to determine hairy vetch aboveground biomass decomposition or N release rates as the model best minimized the standard errors and passed the selection criteria with minimal PRESS and AIC values. However, the modified single exponential model is equally as good as the double exponential for hairy vetch mass and N release rates in the sense that the model can generate valid parameters even for small datasets. To determine the cereal rye decomposition and N release rates, the hyperbolic decay models gave equal weight to the data points, which best estimated the parameters and minimized the residual sum of squares than exponential models. The hyperbolic model had the flexibility to be used with or without asymptote and had a narrow range of rate constants. Hyperbolic models can be used for belowground biomass and can also remove the problem of undue weight to the initial data points that the first-order exponential model gives. Use of best models specific to cereal rye and hairy vetch would give unbiased prediction. We recommend double exponential function with an asymptote for aboveground hairy vetch and the hyperbolic model for cereal rye mass loss and N release. While this modeling study covered hairy vetch and cereal rye decomposition in a humid subtropical region, imposing a different environment might alter the model performance. Specifically, temperature and precipitation may change the rate of decomposition and N release. Investigation of the performance of the best selected models in this study to different legumes and grasses would be the next step in validating the current research results in broader scenarios.