Modeling the Methane Production Kinetics of Anaerobic Co-Digestion of Agricultural Wastes Using Sigmoidal Functions

: The modiﬁed sigmoidal bacteria growth functions (the modiﬁed Gompertz, logistic, and Richards) were used to evaluate the methane production process kinetics of agricultural wastes. The mesophilic anaerobic co-digestion experiments were conducted with various agricultural wastes as feedstocks, including cow manure, corn straw, grape leaves, vines, wine residue, strawberry leaves, and tomato leaves. The results showed that anaerobic co-digestion of cow manure and other agricultural wastes increased the methane yields while it prolonged the lag phase time. Compared with the modiﬁed Gompertz and logistic models, the modiﬁed Richards model obtained higher correlation coefﬁcients and was able to ﬁt experimental data better. The results of this study were expected to determine a suitable model to simulate and study the kinetic process of anaerobic co-digestion with mixed agricultural wastes as feedstocks. Richards model was consistent with that of the modiﬁed Gompertz model. The results suggested that, compared with the modiﬁed Gompertz and Richards models, the ﬁtting result of the CD + VI group obtained from the modiﬁed logistic model was closer to the experimental result.


Introduction
Anaerobic co-digestion of agricultural wastes turns waste into treasure, which not only reduces the pollution of waste to the environment, but also produces economic benefits [1]. Anaerobic co-digestion has attracted many researchers' attention due to its advantages, such as better nutrient balance, higher biogas yield, and higher residual fertilizer value compared with traditional methods of agricultural waste treatment (incineration, landfill, composting, etc.) [2][3][4][5]. Most of the previous studies on anaerobic co-digestion of agricultural wastes focused on a direct analysis of experimental data [6][7][8], instead of modeling the methane production kinetics. In order to understand the impacts of agriculture waste types on anaerobic co-digestion, studying the methane production kinetics of anaerobic co-digestion with different feedstocks is necessary [9].
A bacterial growth curve is used to describe the change in the population of viable bacteria with time (t), and the growth curve represents a growth process with a lag phase just after t = 0 followed by an exponential phase and then by a saturation phase [10]. The cumulative methane yields curves in the early stage (slow methane production), the middle stage (fast methane production, i.e., the peak period), and the later stage (slow methane production and stabilization) of anaerobic co-digestion correspond to the lag, exponential, and saturation phases of the growth curve, respectively, showing a sigmoidal curve [11,12].
Many sigmoidal models were established, such as the Gompertz model, the logistic model, the Richards model, the Stannard model, and the Schnute model [13][14][15][16][17][18]. Among these models, the modified equation of the Stannard model is the same as that of the modified Richards equation. Replacement of a and b (mathematical parameters used in the Schnute equation without biological meanings) in the Schnute equation would result in the modified Richards model [10]. Previous studies showed that the modified Gompertz equation is a common model for methane production by simple organic substrate  1 CD, cow dung; 2 CS, corn straw; 3 GL, grape leaves; 4 VI, vines; 5 WR, wine residue; 6 SL, strawberry leaves; 7 TL, tomato leaves; 8 \, no data; 9 TS, total solid; 10 VS, volatile solid; 11 TOC, total organic carbon; 12 TN, total nitrogen; 13 TP, total phosphorus.
The total solid (TS), volatile solid (VS), soluble sugar, protein, starch, lipid, cellulose, hemicellulose, and lignin were measured using methods previously reported [22]. Total organic carbon (TOC) was measured using the potassium dichromate volumetric method [23]. Total nitrogen (TN) was measured using the indophenol blue colorimetric method after being digested by concentrated sulfuric acid and 30% hydrogen peroxide [24]. Total phosphorus (TP) was measured using a method previously reported [25].
There were nine groups of anaerobic digestion experiments with different agricultural wastes as feedstocks, namely, CD, CD + CS, CD + GL, CD + VI, CD + WR, CD + SL, CD + TL, CD + CS + VI, and CD + CS + WR groups. The compositions and dry weights  Table 2. Among the nine groups, the CD group (with cow dung only) was set as the control group. In addition to the control group, anaerobic co-digestion experiments of two feedstocks were carried out at first. The C/N ratios of CD + CS, CD + VI, CD + SL, and CD + TL groups were 20.61:1, 28.16:1, 20.76:1, and 25.25:1, respectively, which are located in the suitable range of C/N ratio (20-30:1) [26]. In order to keep the same amount of cow dung, the mixture ratio of other groups with two feedstocks (CD + GL and CD + WR) were also set as 2:1. The results for anaerobic co-digestion of both feedstocks showed that the highest A value was obtained in the CD + CS group, followed by the CD + WR group. The combination of CD + CS + WR was then chosen for anaerobic co-digestion to study whether three feedstocks were able to obtain higher biogas yields than the two-feedstock groups. On the other hand, the CD + CS + VI group was chosen to study whether the biogas production can be enhanced and the daily biogas peak can be brought forward compared to the CD + VI group. The C/N ratio of the CD + CS + VI group was 24.38:1, which was suitable for anaerobic digestion. Table 2. The feedstocks and dry weights of the anaerobic co-digestion.

Groups
The Feedstocks and Dry Weights

Experimental Conditions
Before experiments, the feedstocks were mixed according to the dry weights shown in Table 2. The experiments were performed in an Automatic Methane Potential Test System (AMPTS II, Bioprocess Control AB, Lund, Sweden). The TS value of the substrate in the reactor was adjusted to 8% by adding distilled water. Mesophilic (35.0 ± 1.0 • C) experiments were conducted for 60 days in 500 mL fermentation tanks with a working volume of 300 mL. Carbon dioxide was absorbed through carbon dioxide fixation devices, and methane was measured using gas volume measurement devices. Three parallel samples were set for each digestion group. Before analysis, the data of the three parallel samples were averaged and used as the final methane yields. All the methane yields were converted into standard conditions.

Models Description
Assuming that methane production is a function of bacterial growth, i.e., methane production is a composite function, the mathematical parameters in the original function were improved into biological parameters [10], and the modified Gompertz equation [18] can be given by where y represents the cumulative methane yields (mL/g·TS) with respect to time t (d), A is the maximum cumulative methane yields (mL/g·TS), R max is the maximum methane production rate (mL/g·TS/d), λ is the lag phase time (d), and e is a constant equal to 2.71. The logistic model [12] is given by where A, R max , and λ have the same meaning as above. The modified Richards model [10] with the fourth parameter v is given by where A, R max , and λ have the same meaning as above, v is the shape coefficient.

Data Analysis
Nonlinear regression analysis (modified Gompertz, logistic, and Richards) was performed using Origin 9.0 software. The program searched for A, R max , and λ values with the minimum residual sum of squares and their 95% confidence intervals. The curve of methane cumulative yields and R 2 value were obtained during the fitting. R 2 reflected the degree of correlation between variables, and the t-test and F-test were used to conduct statistical tests on the fitting results obtained from different models. According to the analysis of variance (ANOVA) result of the fitting statement, p < 0.05 indicated a significant fitting effect. Figure 1a shows the variation in daily methane yields of anaerobic co-digestion with different feedstocks over time. The daily methane yields in all experimental groups were higher than the control group on the first day of digestion. The control group reached the peak period on the fifth day, while the co-digestion groups attained their peak after the 10th day. From the 26th to 43rd day of the experiment, the daily methane yields in all the experimental groups were generally higher than that in the control group. It was found that anaerobic co-digestion was beneficial to subsequent methane production. Meanwhile, the maximum daily methane yields of experimental groups were also higher than that of the control group. The highest maximum daily methane yield was obtained in the CD + CS group, which was 2.36 times higher than that of the control group.

Methane Yields of Anaerobic Digestion with Different Feedstocks
where y represents the cumulative methane yields (mL/g•TS) with respect to time t (d), A is the maximum cumulative methane yields (mL/g•TS), Rmax is the maximum methane production rate (mL/g•TS/d), λ is the lag phase time (d), and e is a constant equal to 2.71.
The logistic model [12] is given by where A, Rmax, and λ have the same meaning as above. The modified Richards model [10] with the fourth parameter v is given by where A, Rmax, and λ have the same meaning as above, v is the shape coefficient.

Data Analysis
Nonlinear regression analysis (modified Gompertz, logistic, and Richards) was performed using Origin 9.0 software. The program searched for A, Rmax, and λ values with the minimum residual sum of squares and their 95% confidence intervals. The curve of methane cumulative yields and R 2 value were obtained during the fitting. R 2 reflected the degree of correlation between variables, and the t-test and F-test were used to conduct statistical tests on the fitting results obtained from different models. According to the analysis of variance (ANOVA) result of the fitting statement, p < 0.05 indicated a significant fitting effect. Figure 1a shows the variation in daily methane yields of anaerobic co-digestion with different feedstocks over time. The daily methane yields in all experimental groups were higher than the control group on the first day of digestion. The control group reached the peak period on the fifth day, while the co-digestion groups attained their peak after the 10th day. From the 26th to 43rd day of the experiment, the daily methane yields in all the experimental groups were generally higher than that in the control group. It was found that anaerobic co-digestion was beneficial to subsequent methane production. Meanwhile, the maximum daily methane yields of experimental groups were also higher than that of the control group. The highest maximum daily methane yield was obtained in the CD + CS group, which was 2.36 times higher than that of the control group. As can be seen from Figure 1b, the variations of cumulative methane yield over time presented a "stable-rising-stable" trend, showing sigmoidal curves, except for the CD + SL group, CD + TL group, and CD + VI group. The curves of cumulative methane yields of CD + SL and CD + TL groups included multiple "stable-rising" stages. The curve of the CD + VI group presented a phased upward trend. The cumulative methane yields of Energies 2021, 14, 258 5 of 12 CD + CS group were highest, followed by the CD + WR and CD + GL groups, which were 2.37, 2.25, and 2.23 times higher than that of the control group, respectively. The cumulative methane yield of the control group was lowest. According to the daily and cumulative methane yields, it was found that anaerobic co-digestion increased the cumulative methane yields by increasing the maximum daily methane yields and expanding the period with high yields.

Modeling Methane Production Kinetics Using Modified Gompertz Model
The cumulative methane yields obtained from the anaerobic co-digestion experiments were fitted using the modified Gompertz, logistic, and Richards models ( Figure 2). The comparison between fitting results obtained from the modified Gompertz model and experimental results is shown in Figure 2a. As can be seen from the data in Table 3, the R 2 values were all greater than 0.9, which showed the good fitting results of the modified Gompertz model in expressing the accumulative process of methane yields. The highest R 2 values were obtained in the CD group. The R 2 value of the CD + TL group was the lowest.
As can be seen from Figure 1b, the variations of cumulative methane yield over time presented a "stable-rising-stable" trend, showing sigmoidal curves, except for the CD + SL group, CD + TL group, and CD + VI group. The curves of cumulative methane yields of CD + SL and CD + TL groups included multiple "stable-rising" stages. The curve of the CD + VI group presented a phased upward trend. The cumulative methane yields of CD + CS group were highest, followed by the CD + WR and CD + GL groups, which were 2.37, 2.25, and 2.23 times higher than that of the control group, respectively. The cumulative methane yield of the control group was lowest. According to the daily and cumulative methane yields, it was found that anaerobic co-digestion increased the cumulative methane yields by increasing the maximum daily methane yields and expanding the period with high yields.

Modeling Methane Production Kinetics Using Modified Gompertz Model
The cumulative methane yields obtained from the anaerobic co-digestion experiments were fitted using the modified Gompertz, logistic, and Richards models ( Figure 2). The comparison between fitting results obtained from the modified Gompertz model and experimental results is shown in Figure 2a. As can be seen from the data in Table 3, the R 2 values were all greater than 0.9, which showed the good fitting results of the modified Gompertz model in expressing the accumulative process of methane yields. The highest R 2 values were obtained in the CD group. The R 2 value of the CD + TL group was the lowest.  As shown in Table 3, for the CD group and CD + TL group, the A values of fitting results were lower than those of experimental results. For other groups, the fitted A values were higher than those obtained from experiments. The A values with high experimental results obtained high fitting results, except for the CD + VI group. The results showed that A value obtained from the model of the CD + VI group was quite different from the experimental result. This indicated that the cumulative methane yield curve of the CD + VI group was not a sigmoidal curve, and the cumulative methane yield of this group did not conform to the three-stage hypothesis of the modified Gompertz function.
The R max values of the experimental results were higher than those of the fitting results in all groups ( Table 3). The highest and lowest fitting R max values were obtained in the CD + GL group and CD + SL group, respectively. However, the highest and lowest experimental R max values were obtained in the CD + CS group and CD group, respectively. Therefore, the fitting results were not in line with the experimental results.
The lag phase times (λ) of experimental groups were 99.6-189% longer than those of the control group with the exception of the fitting result obtained in the CD + SL group. The λ values of fitting results in these groups were more than 10 days. In a mathematical sense, λ represents the point at which it intersects the x-axis. Biologically speaking, λ is the time of the slow methane production stage. Therefore, the fitting results were consistent with the experimental results. The modified Gompertz model did not provide a very accurate fit to the experimental data of the CD + SL group, even though the R 2 value was 0.9503. The lag phase time (λ) of this group was 3.86 days (Table 3), which was 23.1% shorter than that of the control group. The daily methane yields of the CD + SL group were low from the third to the 10th day of the experiment (Figure 1a). Therefore, the reasonable λ value of the CD + SL group was more than 10 days. The inadequate fit of the modified Gompertz equation was understood for the particular features of the original growth function and the curve of this group was not sigmoidal. As the point of inflection was fixed for the modified Gompertz function, the remainder of the curve was modeled in accordance with the location point, as well as the asymmetry for the Gompertz [12].

Modeling Methane Production Kinetics Using Modified Logistic Model
The cumulative methane yields of anaerobic co-digestion fitted by the modified logistic equation and experimental results are shown in Figure 2b. As can be seen from the data in Table 4, R 2 values obtained from modified logistic model were higher than 0.9, indicating the good fit of the logistic model in expressing the accumulation process of methane yields. R 2 values fitted by the modified logistic model were higher than those fitted by the modified Gompertz model, except for the CD group and CD + CS + WR group. For the CD + CS group, both the fitting result and the experimental result of the A value were highest (Table 4), which was inconsistent with the modified Gompertz model. As can be seen from the data in Tables 3 and 4, the A values obtained from the modified logistic model were lower than those obtained from the modified Gompertz model. The A values obtained from the modified logistic model was closer to the A values of the experiment compared to the modified Gompertz model, with the exception of the CD group and CD + TL group.
The R max values of the experimental results were higher than those of the fitting results (Table 3). However, the highest or lowest fitting R max values and experimental R max values were not obtained in the same group. The reason was that the R max values of experimental results were determined by calculating the value of maximum methane production rates. The program estimated the R max values by searching for the steepest ascent of the curve [10], i.e., the R max values were determined from the slope of the line during exponential gas production incorporating multiple data points [12]. Therefore, the fitting R max values were not completely consistent with the experimental R max values.
The λ values of experimental groups obtained from the logistic model were higher than those of the control group (Table 4). The λ values of experimental groups fitted by the modified logistic model were higher than those fitted by the modified Gompertz model. The differences in results between these two models were interpreted as being due to the different R max values estimated by the program, leading to different λ values, as the λ values were determined by intersecting the tangent to the point of inflection with the x-axis [10,27].

Modeling Methane Production Kinetics Using Modified Richards Model
The modified Richards equation was used to simulate the cumulative methane yields, and the results are shown in Figure 2c. The average R 2 values obtained from the modified Richards model were higher than those obtained from the modified Gompertz and logistic models. The results showed that, compared with the modified Gompertz and logistic models, better fitting results were obtained in the modified Richards model.
The highest A values of experimental results did not correspond to the highest A values of fitting results (Table 5). For the relationship between A values of fitting results and experimental results, the conclusion of the modified Richards model was consistent with that of the modified Gompertz model. The results suggested that, compared with the modified Gompertz and Richards models, the fitting result of the CD + VI group obtained from the modified logistic model was closer to the experimental result. The R max values of experimental results were higher than those of fitting results, with the exception of the CD + WR group. The fitting R max value of the CD + WR group was highest. The fitting A value of the CD + VI group was highest. The R max values did not correspond to the A values, which were consistent with results of the modified Gompertz and logistic models. It was understood that R max values represent the maximum daily methane yields when the curves were the standard sigmoidal curves. Some of the curves in this study were not standard sigmoidal curves, such as the curves of the CD + VI group, CD + SL group, and CD + TL group.
The λ values of experimental groups fitted by the modified Richards model were higher than those of the control group, which was consistent with the results obtained from the modified logistic model. The results showed that anaerobic co-fermentation with different agricultural wastes prolonged the start-up period of fermentation. For the λ value of the CD + SL group, the modified Gompertz and logistic models failed to accurately represent the experimental data, and the modified Richards model provided a relatively good fit to the experimental data. The reason was that the curve of this group was not sigmoidal, and different models (the points of inflection were fixed for the modified Gompertz and logistic model, while the point of inflection was flexible for the modified Richards model) provided different results when estimated by the algorithm [12]. Therefore, compared with the modified Gompertz and logistic models, the modified Richards model obtained a more accurate result in the CD + SL group, which indicated that the Richards model was more adaptable.

Discussion
According to the results shown in Figure 1, it was found that anaerobic co-digestion obtained higher maximum daily and cumulative methane yields than the group with cow dung only. On one hand, the C/N ratios in most groups were suitable for anaerobic digestion compared with the control group. On the other hand, anaerobic co-digestion balanced the nutrients better and was beneficial to the methane production. The anaerobic co-digestion reached its peak period later. It was speculated that the complexity of mixed substrates affected the start-up time of anaerobic co-digestion.
As shown in Tables 3-5, the highest R 2 value was obtained from the modified Gompertz model in the control group, which was consistent with previous studies [10,28,29]. The R 2 values obtained from the modified Richards model were higher than 0.99, with the exception of the CD + SL group and CD + TL group. It was speculated that the CD + SL group and CD + TL group basically did not complete methane production. The basic completion of anaerobic digestion is based on the fact that the daily methane yields are less than 1% of the cumulative methane yields [29]. The average daily methane yields of the CD + SL group and CD + TL group were 0.42 mL/g·TS and 1.67 mL/g·TS, while the maximum daily Energies 2021, 14, 258 9 of 12 methane yields of these two groups were 27.19 mL/g·TS and 35.83 mL/g·TS, which were 1.57% and 4.67% of the maximum daily methane yields, respectively. Therefore, anaerobic digestion of the CD + SL group and CD + TL group was not completed. Therefore, the curves of these two groups did not show complete sigmoidal curves, resulting in lower R 2 values than other groups.
The average R 2 values of the modified Richards model were higher than those of the modified Gompertz and logistic models. The correlation coefficient (R 2 ) is a statistical index used to reflect the degree of correlation between variables. In order to verify the significance of these three models in the anaerobic co-digestion with agricultural wastes as feedstocks, a t-test and F-test were conducted. Statistical analysis of the fitting report showed that the parameters in the modified Gompertz, logistic, and Richards models were suitable for estimating the experimental results, except for the R max and v values of the CD group and CD + CS + WR group and the v value of the CD + VI group obtained from the modified Richards model. At the 0.05 level, the fitting function was significantly better than the function of experimental data.
Although the R 2 values of the three models were relatively high and the fitting curves ( Figure 2) of the three models were consistent with the experimental results on a whole, the fitting results of A values and R max values of the three models were different. Figure 3 shows the difference among the A values fitted by the three models and the experimental results. The results suggested that, when the A values of experimental results were high, the A values of fitting results were also high with the exception of the CD + VI group. For the CD + VI group, the A values obtained from the three models were higher than the A values of the experiment results. This is understandable because the curves of the CD + VI group fitted by the three models were not sigmoidal. The fitted curves showed an upward trend, i.e., the fitted A values appeared after the 60th day, which were higher than the experimental A values.
Energies 2021, 14, x FOR PEER REVIEW 9 of 12 group and CD + TL group basically did not complete methane production. The basic completion of anaerobic digestion is based on the fact that the daily methane yields are less than 1% of the cumulative methane yields [29]. The average daily methane yields of the CD + SL group and CD + TL group were 0.42 mL/g•TS and 1.67 mL/g•TS, while the maximum daily methane yields of these two groups were 27.19 mL/g•TS and 35.83 mL/g•TS, which were 1.57% and 4.67% of the maximum daily methane yields, respectively. Therefore, anaerobic digestion of the CD + SL group and CD + TL group was not completed. Therefore, the curves of these two groups did not show complete sigmoidal curves, resulting in lower R 2 values than other groups. The average R 2 values of the modified Richards model were higher than those of the modified Gompertz and logistic models. The correlation coefficient (R 2 ) is a statistical index used to reflect the degree of correlation between variables. In order to verify the significance of these three models in the anaerobic co-digestion with agricultural wastes as feedstocks, a t-test and F-test were conducted. Statistical analysis of the fitting report showed that the parameters in the modified Gompertz, logistic, and Richards models were suitable for estimating the experimental results, except for the Rmax and v values of the CD group and CD + CS + WR group and the v value of the CD + VI group obtained from the modified Richards model. At the 0.05 level, the fitting function was significantly better than the function of experimental data.
Although the R 2 values of the three models were relatively high and the fitting curves ( Figure 2) of the three models were consistent with the experimental results on a whole, the fitting results of A values and Rmax values of the three models were different. Figure 3 shows the difference among the A values fitted by the three models and the experimental results. The results suggested that, when the A values of experimental results were high, the A values of fitting results were also high with the exception of the CD + VI group. For the CD + VI group, the A values obtained from the three models were higher than the A values of the experiment results. This is understandable because the curves of the CD + VI group fitted by the three models were not sigmoidal. The fitted curves showed an upward trend, i.e., the fitted A values appeared after the 60th day, which were higher than the experimental A values. The fitting results of the three models showed that there was no clear relationship between the A values and the R 2 values (Tables 3-5 and Figure 3), which was consistent with previous research [11]. The high R 2 values indicated that the degree of correlation between variables was high. The A values were related to the type of feedstock, proportion of different feedstocks, and conditions of anaerobic fermentation. Therefore, the conclusion agreed with practice.  The fitting results of the three models showed that there was no clear relationship between the A values and the R 2 values (Tables 3-5 and Figure 3), which was consistent with previous research [11]. The high R 2 values indicated that the degree of correlation between variables was high. The A values were related to the type of feedstock, proportion of different feedstocks, and conditions of anaerobic fermentation. Therefore, the conclusion agreed with practice.
The results of ANOVA showed that there was no significant difference for A values between the fitting values obtained from the three models and the experimental values with the exception of the CD + VI group (Figure 3). The multiple daily methane yield peaks in the CD + VI group (Figure 1a), which might be caused by the different properties of VI, resulted in the inaccuracy of fitting. Therefore, the modified Gompertz, logistic, and Richards models were all able to represent the experimental A values in most of the groups except for the CD + VI group. The A values of the CD + VI group fitted by the logistic and Richards models were closer to the A values of the experimental results than the Gompertz model. Therefore, the A values of the CD + VI group obtained from the modified logistic and Richards models were more representative of the experimental A values than those obtained from the modified Gompertz model. Figure 4 shows the differences between the experimental values and the fitting values of R max . The results showed that the fitted R max values of the three models were lower than the experimental R max values except for the CD + WR group. Lower fitted R max values than the experimental R max values were also reported in previous anaerobic co-digestion studies [12,30]. The average R max values fitted by the modified logistic and Richards models were closer the experimental R max values than those fitted by the modified Gompertz model. Therefore, for the sake of R max , the modified logistic and Richards models represented the experimental results better than the modified Gompertz models.
The results of ANOVA showed that there was no significant difference for A values between the fitting values obtained from the three models and the experimental values with the exception of the CD + VI group (Figure 3). The multiple daily methane yield peaks in the CD + VI group (Figure 1a), which might be caused by the different properties of VI, resulted in the inaccuracy of fitting. Therefore, the modified Gompertz, logistic, and Richards models were all able to represent the experimental A values in most of the groups except for the CD + VI group. The A values of the CD + VI group fitted by the logistic and Richards models were closer to the A values of the experimental results than the Gompertz model. Therefore, the A values of the CD + VI group obtained from the modified logistic and Richards models were more representative of the experimental A values than those obtained from the modified Gompertz model. Figure 4 shows the differences between the experimental values and the fitting values of Rmax. The results showed that the fitted Rmax values of the three models were lower than the experimental Rmax values except for the CD + WR group. Lower fitted Rmax values than the experimental Rmax values were also reported in previous anaerobic co-digestion studies [12,30]. The average Rmax values fitted by the modified logistic and Richards models were closer the experimental Rmax values than those fitted by the modified Gompertz model. Therefore, for the sake of Rmax, the modified logistic and Richards models represented the experimental results better than the modified Gompertz models. The modified Gompertz, logistic, and Richards models showed similar results for the λ values with the exception of the CD + SL group, i.e., the λ values of the experimental groups were higher than those of the control group. The study by Ware et al. suggested that low λ values demonstrated the simple nature and high biodegradability of the feedstocks [12]. Because the anaerobic co-digestion experiment of agricultural wastes contained a higher content of lignocellulose than the experiment of the control group (with cow manure only), the hydrolysis rate was limited [31,32]. Thus, the digestion times were affected. The results showed that, although the cumulative methane yields were increased by anaerobic co-digestion, the start-up period was extended. For the CD + SL group, compared with the results obtained from the modified Gompertz and logistic models, λ values fitted by the modified Richards model better represented experimental results (Tables 3-5). It was speculated that the modified Richards model is a four-parameter model, and the fourth parameter v was introduced to make the shape of the curve flexible. The fitting result of the four-parameter model was more reasonable than that of the three-parameter models (modified Gompertz and logistic models). Moreover, comparing the results of R 2 , A, Rmax, and λ, it was found that all the parameters should be taken into account when evaluating the accuracy of models.  The modified Gompertz, logistic, and Richards models showed similar results for the λ values with the exception of the CD + SL group, i.e., the λ values of the experimental groups were higher than those of the control group. The study by Ware et al. suggested that low λ values demonstrated the simple nature and high biodegradability of the feedstocks [12]. Because the anaerobic co-digestion experiment of agricultural wastes contained a higher content of lignocellulose than the experiment of the control group (with cow manure only), the hydrolysis rate was limited [31,32]. Thus, the digestion times were affected. The results showed that, although the cumulative methane yields were increased by anaerobic co-digestion, the start-up period was extended. For the CD + SL group, compared with the results obtained from the modified Gompertz and logistic models, λ values fitted by the modified Richards model better represented experimental results (Tables 3-5). It was speculated that the modified Richards model is a four-parameter model, and the fourth parameter v was introduced to make the shape of the curve flexible. The fitting result of the four-parameter model was more reasonable than that of the three-parameter models (modified Gompertz and logistic models). Moreover, comparing the results of R 2 , A, R max , and λ, it was found that all the parameters should be taken into account when evaluating the accuracy of models.
In general, the fitting results of the modified Gompertz model were better in expressing the methane production process of the control group. The modified logistic model was able to replicate the data of experimental groups. The modified Richards model performed better than the other two models in fitting the experimental results of anaerobic co-digestion with different feedstocks. The results of this study were consistent with a previous study [12]. The three-parameter sigmoidal model was insufficient in accurately describing methane yields from complex feedstocks. Due to the nature of biological models, there is no one-size-fits all dynamic model to fit the cumulative methane yields with different feedstocks. In order to determine the most suitable model, three-parameter and four-parameter models should be applied simultaneously [10,33].

Conclusions
This study showed that, compared with the anaerobic digestion of cow dung only, the anaerobic co-digestion of agricultural wastes (two or three feedstocks) increased the cumulative methane yields, while prolonged the lag phase time. According to the R 2 values and the differences between the experimental and fitting values, the better fitting results of the experimental groups were obtained from the modified Richards model compared with those obtained from the modified Gompertz and logistic models. In order to determine the most suitable model, three-parameter and four-parameter models should be applied simultaneously.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data are presented in this article. Data sharing is not applicable to this article.