Effects of Organic Energy Crop Rotations and Fertilisation with the Liquid Digestate Phase on Organic Carbon in the Topsoil

: Combining organic farming and biogas production from agricultural feedstocks has been suggested as a way of achieving carbon (C) neutrality in Europe. However, as the long-term effects of C removal for methane production on soil organic carbon (SOC) are unclear, organic farmers in particular have questioned whether farm biogas production will have a positive effect on soil fertility. Eight years of data from an organic long-term ﬁeld trial involving digestate fertilisation and various crop rotations (CRs) with differing proportions of clover-grass leys were used to calculate C inputs based on the CANDY model, and these modelled changes compared with measured changes in SOC content (SOCc) over the same period. Measured SOCc increased by nearly 20% over the eight years. Digestate fertilisation signiﬁcantly increased SOCc. Fertilised plots with the highest proportion of clover-grass in the CR had the highest SOCc. The C inputs from clover-grass leys, even if they only made up 25% of the CR, were high enough to increase SOCc, even with the removal of all aboveground biomass and without fertilisation. Our results show that biogas production based on clover-grass leys could be an important part of sustainable farming, improving or maintaining SOCc and improving nutrient ﬂows, particularly in organic farming, while simultaneously providing renewable energy.


Introduction
Organic farmers in particular are reliant on organic inputs as nutrient sources, as the use of synthetic chemical fertilisers is not permitted in organic farming (see, e.g., The Council of the European Union [1]). Organic fertilisers, however, are usually available in lower amounts than synthetic fertilisers. The nutrients in organic fertilisers often have to be mineralised in soils first before they can be used by crops [2]. As a result, organic farmers are therefore also reliant on soil organic matter acting as a reservoir for nutrients necessary for crop growth [3]. It is, hence, particularly important in organic agriculture that soil organic matter contents are, at the very least, maintained. Preservation of soil organic matter also has other benefits such as improved soil structure and therefore increased water infiltration and storage, and lower erosion risk. Improving soil quality is hence an important part of climate change adaptation in agriculture. Carbon (C) sequestration consequently contributes to food security by improving soil quality [4].
C sequestration in soils is itself essential for climate change mitigation, as is the transition to renewable energies. A scenario for achieving C neutrality in Europe developed by Aubert et al. [5] recommends reducing livestock numbers and instead using grassland biomass for biogas production. As the proportion of farms without livestock is increasing in organic farming in Germany [6], a similar approach could also appeal to organic farmers. Harvesting the biomass from forages such as clover-grass leys for biogas production instead of mulching provides an additional source of income and has been shown to be studies use digestates based on liquid animal manures and/or maize silage, whereas the digestate used in our field trial is produced by an organic farm and is based on clover-grass and grass silage, and therefore more accurately represents the types of digestates that could be used in future sustainable farming scenarios [5].

Experimental Site
Our field trial ( Figure 1) was situated at the Technical University of Munich's (TUM) agricultural research station in Viehhausen. Viehhausen is located approximately 8 km west of Freising in southern Germany, 490 m above sea level in the Tertiary hill country, an undulating landscape developed in unconsolidated Tertiary sediments and overlain by a thin loess cover. Using the US soil taxonomy [36], the soil at the experimental site is categorised as a Hapludalf derived from loess with silty loam texture down to at least 1 m. Using the World Reference Base [37], the soil is a Haplic Luvisol (Manganiferric, Siltic). Average clay content was 25%. The experiment was situated on a slope facing north-east with a gradient of about 9%. Average annual temperature and precipitation for the period 2010-2017 in Viehhausen were 9.0 • C and 799 mm a −1 , respectively. The field trial began in late 2004 and has been running in its present form since 2007.
Agronomy 2021, 11, x FOR PEER REVIEW 3 of 18 work that substantially modify short-term responses in the long run. In addition, many of these studies use digestates based on liquid animal manures and/or maize silage, whereas the digestate used in our field trial is produced by an organic farm and is based on clovergrass and grass silage, and therefore more accurately represents the types of digestates that could be used in future sustainable farming scenarios [5].

Experimental Site
Our field trial ( Figure 1) was situated at the Technical University of Munich's (TUM) agricultural research station in Viehhausen. Viehhausen is located approximately 8 km west of Freising in southern Germany, 490 m above sea level in the Tertiary hill country, an undulating landscape developed in unconsolidated Tertiary sediments and overlain by a thin loess cover. Using the US soil taxonomy [36], the soil at the experimental site is categorised as a Hapludalf derived from loess with silty loam texture down to at least 1 m. Using the World Reference Base [37], the soil is a Haplic Luvisol (Manganiferric, Siltic). Average clay content was 25%. The experiment was situated on a slope facing north-east with a gradient of about 9%. Average annual temperature and precipitation for the period 2010-2017 in Viehhausen were 9.0 °C and 799 mm a −1 , respectively. The field trial began in late 2004 and has been running in its present form since 2007.

Digestate and Fertilisation
The main factor studied in the field trial was fertilisation using biogas digestate, which was assigned in two levels to plots within the factor CR; there were ten different

Digestate and Fertilisation
The main factor studied in the field trial was fertilisation using biogas digestate, which was assigned in two levels to plots within the factor CR; there were ten different CRs and the 32 plots of each CR were divided into 16 fertilised and 16 unfertilised plots ( Figure 1). Each unfertilised plot was neighboured laterally by a fertilised plot on one side across the slope and on one side along the slope. There was no difference in treatment between fertilised and unfertilised plots except that fertilised plots were fertilised using biogas digestate. Each CR received a different amount of digestate, depending on the theoretical amount of digestate that could be produced from the estimated crop biomass of clovergrass, rye (wholecrop) (Secale cereale L.), silage maize (Zea mays L.), triticale (wholecrop) (x Triticosecale Wittmack) and sunflower (wholecrop) (Helianthus annuus L.) grown in the fertilised plots of each CR. The total amount of digestate was then divided among the crops in the CR according to estimated nutrient requirements. The amount of C added via digestate to each plot therefore ranged from 0 to 3 t ha −1 a −1 . The digestate was applied using a slurry tanker fitted with trailing hoses. The digestate was produced by a local organic farmer from a feedstock mix of, on average, 61% silage from clover-grass ley/grass ley/grassland biomass, 30% solid cattle manure, 5% silage maize and the remainder cereal grains. From 2010 onwards, the digestate was separated into liquid and solid phases and the liquid phase used in this field trial (for the sake of simplicity, we refer only to digestate). The average dry matter content of the digestate after separation of the solids was 7.9%, with a total C content of 39% (

Experimental Design
The factor CR comprised 10 different four-year crop rotations ( Table 2). The CRs were laid out in a replicated control design. CR1a was located at the start of the field, repeated in the middle (CR1b) and at the opposite end of the field (CR1c) in order to capture any trends in yield potential across the site. Each crop of every CR was cultivated each year; the CRs were therefore laid out in columns running down the field, divided into four blocks. The crops moved up a block each year. Within each CR there were 8 replications in each block (4 fertilised and 4 unfertilised plots). Plots measured 6 m × 12 m (plots on both outer edges of the trial measured 9 m × 12 m), with alleys of 3 m between each row of plots and 9 m between each block. The field trial therefore had a total size of nearly 4 ha.
The first year in each CR was a lucerne-clover-grass ley or clover-grass ley. The clovergrass leys were mixtures of red clover (both diploid and tetraploid, Trifolium pratense L.) and white clover (Trifolium repens L.) with various grasses, such as orchard grass (Dactylis glomerata L.), meadow fescue (Festuca pratensis Huds.) or Italian ryegrass (Lolium multiflorum Lam.). The lucerne-clover-grass leys also included lucerne (Medicago sativa L.). For the sake of simplicity, we refer only to clover-grass. The ley year was followed by a year of winter wheat (Triticum aestivum L.); the CRs differed in the third and fourth years. All aboveground biomass was harvested and removed (including wheat straw), with the following exceptions: field bean (Vicia faba L.) and soybean (Glycine max (L.) Merr.) straw in CRs 3 and 4; year 3 cover crops in CRs 3, 4, 7 and 8 (mulched); white clover in CRs 5 and Agronomy 2021, 11, 1393 5 of 18 6 (mulched). Leys were cut 3-5 times per year. Conventional tillage (ploughing to a depth of approximately 22 cm) was used. As ploughing depended on which crop was grown, plots were not ploughed every year.

CR1c
Lucerneclover-grass Wheat Biomass and grain yields (Table 3) were determined each year using a plot harvester and plot combine, respectively. Straw and cover crop biomass yields were also measured, but not in all years. The C content of soil and crop grain and biomass samples was analysed using a C and N analyser (Vario Max, Elementar Analysensysteme GmbH, Dumas combustion method). The mean crop C content was between 43 and 45%, with the exception of soybean, where the mean C content was 52%. For the 2017 soil samples, if the Agronomy 2021, 11, 1393 6 of 18 pH indicated the presence of carbonates (pH > 6.8), carbonate-C content was additionally determined after treatment in a muffle furnace at 550 • C [38]. The pH of the soil samples was measured in a 0.01 M CaCl 2 solution using a pH meter. Soil bulk density was measured in August 2014 in block 2 and CRs 8, 10 and 1c, using intact soil cores to a depth of 5 cm. Mean bulk density was 1.4 g cm −3 [39], which was used for the calculation of the size of SOC stocks.

Modelling
C input (kg ha −1 ) for the period 2010 to 2017 was calculated using Equation (1).
C RHR = C input from root and harvest residues (kg ha −1 ) C S = C input from seed (kg ha −1 ) C SF = C input from straw fertilisation (kg ha −1 ) C GM = C input from green manure (kg ha −1 ) C D = C input from digestate (kg ha −1 ). As the C content of seed was not measured, values for C S were taken from the REPRO model [40]. C SF , C GM and C D were calculated based on field trial data for yield, grain and biomass C content, and digestate amount and composition data, for the period 2010-2017. C RHR was calculated using yield data and Equation (2). Equation (2) was taken from the CANDY model [34].
C RHR = C input from root and harvest residues (dt ha −1 ) The crop parameters K RHR and F RHR (Table 4) were adapted from the CANDY model [34] as follows. The averages of the individual values for lucerne, clover and grass were used as the K RHR values for lucerne-clover-grass and clover-grass. The F RHR value for lucerne was used as the F RHR value for lucerne-clover-grass. The crop parameters for soybean, sorghum and triticale were assumed to be the same as the crop parameters for field bean, silage maize and rye, respectively. The crop parameters for wholecrop sunflower were derived from the parameters for grain sunflower, using the ratio of the parameters for grain maize and silage maize. As C RHR for undersown clover-grass was not available in the CANDY model, the value from REPRO was used. Table 4. Crop-specific C input parameters K RHR , F RHR and synthesis coefficient S used to calculate C RHR and C input to soil organic carbon (SOC), adapted from the CANDY model [34]. The C input to SOC (the C that remains after respiration and therefore contributes to SOC stock) was calculated using C RHR , C SF , C GM and C D and the synthesis coefficients S (Table 4), according to Equation (3). Equation (3) and the values for S were taken from the CANDY model [34].
It was assumed that C S was used by the growing crop and therefore did not contribute to SOC stocks. The S for the digestate used in our field trial was assumed to be the same as the S for animal-based organic fertilisers (e.g., solid or liquid manure), and the S for bean straw was assumed to be the same as S for cereal straw.

Statistical Analysis
Statistical analyses were carried out using R version 3.5.2 [41]. To analyse the differences in SOCc between years and fertilised and unfertilised plots, a linear mixed model using restricted maximum likelihood estimation criterion (REML) with the fixed effects year, fertilisation and block was fitted. Intercepts were allowed to vary for the random effect CR field position along the x-axis (n = 12). To analyse the differences in SOCc between years and CRs, a linear mixed model using REML with the fixed effects year, CR and block was fitted, and field position of each plot along the x-axis was used as the random factor. A type III ANOVA with Satterthwaite degrees of freedom approximation was calculated, followed by Tukey post hoc tests of the estimated marginal means. The R packages lmer4 and emmeans were used for these analyses. A two-way factorial ANOVA was used to analyse the differences in SOCc between CRs 1a, 1b and 1c, with the factors year and fertilisation. Linear regression was used to analyse the relationships between digestate fertilisation, CR components and SOCc.

Results
SOCc increased significantly from 1.10% in 2010 to 1.31% in 2017 (t(378) = 30.94, p < 0.001). This increase was significant in both the fertilised and the unfertilised plots (Table 5; for illustration see Figure 2). The increase was larger in the fertilised plots than in the unfertilised plots (23% versus 15% of the SOCc in 2010).
In 2010, six years after the start of the experiment, SOCc in the fertilised plots was 4% higher than in the unfertilised plots and, in 2017, 11% higher, indicating that the effect had increased over time. . This increase was significant in both the fertilised and the unfertilised plots (Table  5; for illustration see Figure 2). The increase was larger in the fertilised plots than in the unfertilised plots (23% versus 15% of the SOCc in 2010). In 2010, six years after the start of the experiment, SOCc in the fertilised plots was 4% higher than in the unfertilised plots and, in 2017, 11% higher, indicating that the effect had increased over time. An increase occurred in nearly all plots (all data points above the 1:1 line in Figure 2) but the scatter was considerable when looking at individual (unreplicated) plots, with random variation in both the 2010 and in the 2017 measurements.
In 2010, SOCc was highest in CR6 with 50% maize undersown with white clover in the rotation, and lowest in CR2 and CR3, both with two years of wheat in the rotation. SOCc increased significantly (p < 0.001) in all CRs between 2010 and 2017 ( Figure 3). In 2017 in both the fertilised and unfertilised plots, SOCc was highest in CR10, with 75% clover-grass ley in the rotation, and lowest in CR3. Digestate fertilisation was lowest in CR3. CR10 also had the biggest increase in SOCc between 2010 and 2017 (26%), whereas in CR4, with two years of wheat and one year of soybean in the rotation, SOCc increased by only 13%. In the fertilised plots the increase was lowest in CR6; however, these plots An increase occurred in nearly all plots (all data points above the 1:1 line in Figure 2) but the scatter was considerable when looking at individual (unreplicated) plots, with random variation in both the 2010 and in the 2017 measurements.
In 2010, SOCc was highest in CR6 with 50% maize undersown with white clover in the rotation, and lowest in CR2 and CR3, both with two years of wheat in the rotation. SOCc increased significantly (p < 0.001) in all CRs between 2010 and 2017 ( Figure 3). In 2017 in both the fertilised and unfertilised plots, SOCc was highest in CR10, with 75% clover-grass ley in the rotation, and lowest in CR3. Digestate fertilisation was lowest in CR3. CR10 also had the biggest increase in SOCc between 2010 and 2017 (26%), whereas in CR4, with two years of wheat and one year of soybean in the rotation, SOCc increased by only 13%. In the fertilised plots the increase was lowest in CR6; however, these plots did have the highest mean SOCc in 2010. In CR9 SOCc was 14% higher in 2017 in the fertilised plots than in the unfertilised plots, whereas in CR4 SOCc in the fertilised plots was only 6.5% higher. CR9 received the largest amount of digestate in the experiment, whereas CR4 received the second-lowest amount. Despite this, the difference between the fertilised and unfertilised plots changed (increased) the most in CR4. In comparison, the relative difference between the fertilised and unfertilised plots stayed almost constant in CR6.
Despite identical treatment, SOCc was significantly lower in CR1a than in CR1b or CR1c in both 2010 and 2017, for both the fertilised and unfertilised plots, p < 0.01.
Agronomy 2021, 11, 1393 9 of 18 did have the highest mean SOCc in 2010. In CR9 SOCc was 14% higher in 2017 in the fertilised plots than in the unfertilised plots, whereas in CR4 SOCc in the fertilised plots was only 6.5% higher. CR9 received the largest amount of digestate in the experiment, whereas CR4 received the second-lowest amount. Despite this, the difference between the fertilised and unfertilised plots changed (increased) the most in CR4. In comparison, the relative difference between the fertilised and unfertilised plots stayed almost constant in CR6.  Table 6 for significant differences between CRs.   Table 6 for significant differences between CRs. There was a significant correlation between the proportion of clover-grass leys on a plot and the SOCc (Figure 4). SOCc in 2017 was higher, for both the fertilised and unfertilised plots, in those plots that had had more years of clover-grass leys. Correspondingly, CR10 (75% clover-grass leys in the CR) and CR9 (50% clover-grass leys in the CR) had the highest SOCc in 2017 (Figure 4). The correlation was slightly stronger for the fertilised plots, however the difference between the slopes of the regressions for the fertilised and unfertilised plots was not significant. This indicates that the effect of the clover-grass leys on SOCc in 2017 was independent of the effect of digestate. Using a combined slope for the fertilised and unfertilised plots and keeping the fertiliser effect constant yielded the following equation that, although it includes only two parameters, explained 44% of the variation: F = Fertiliser total C amount (kg ha −1 ) L = Ley proportion (%) R 2 = 0.44, n = 379. Both effects were very highly significant (p < 0.001).
Despite identical treatment, SOCc was significantly lower in CR1a than in CR1b or CR1c in both 2010 and 2017, for both the fertilised and unfertilised plots, p < 0.01.
There was a significant correlation between the proportion of clover-grass leys on a plot and the SOCc (Figure 4). SOCc in 2017 was higher, for both the fertilised and unfertilised plots, in those plots that had had more years of clover-grass leys. Correspondingly, CR10 (75% clover-grass leys in the CR) and CR9 (50% clover-grass leys in the CR) had the highest SOCc in 2017 (Figure 4). The correlation was slightly stronger for the fertilised plots, however the difference between the slopes of the regressions for the fertilised and unfertilised plots was not significant. This indicates that the effect of the clover-grass leys on SOCc in 2017 was independent of the effect of digestate. Using a combined slope for the fertilised and unfertilised plots and keeping the fertiliser effect constant yielded the following equation that, although it includes only two parameters, explained 44% of the variation: SOCc = 1.13 + 1.8 × 10 −5 × F + 0.003 × L,  In the CRs with only one year of clover-grass there was a significant positive correlation between SOCc in 2017 and the proportion of row crops with undersown crops ( Figure   Figure 4. Correlation between the proportion of clover-grass leys on a plot in the period 2005-2017, and the SOCc (%) of that plot in 2017. Unfertilised plots are represented by grey triangles, fertilised plots by black circles. The dashed line is the combined regression line (y = 1.19 + 0.004x, R 2 = 0.13, p < 0.001, n = 379). The grey solid line is the regression line for the unfertilised plots (y = 1.13 + 0.004x), the black solid line is the regression line for the fertilised plots (y = 1.26 + 0.004x), assuming no interaction between the proportion of clover-grass and fertilisation. Jitter was added to the proportion of clover-grass leys to improve visibility of symbols.
In the CRs with only one year of clover-grass there was a significant positive correlation between SOCc in 2017 and the proportion of row crops with undersown crops ( Figure 5). Note that, for the sake of simplicity, we refer to the white clover in CRs 5 and 6 as an undersown crop, although the maize was sown into the clover cover crop. This correlation was significant for both the fertilised and unfertilised plots (p < 0.05) but again, the slopes of the regressions for the fertilised and unfertilised plots were not significantly different (p = 0.19). A combined slope could therefore be assumed, which resulted in the following equation: U = Proportion undersown row crops (%) R 2 = 0.38, n = 316. SOCc = 1.11 + 1.7 × 10 −5 × F + 0.004 × L + 0.002 × U, R 2 = 0.46, n = 379. All effects were very highly significant (p < 0.001). R 2 indicated that these three influences explained 46% of the variation in SOCc, while 54% can be attributed to random variation and other influences (e.g., crop type, interactions). Mean total C input, as calculated using Equation (1), was 1.98 t ha −1 a −1 in the unfertilised plots and 3.17 t ha −1 a −1 in the fertilised plots. Of this, 0.89 t ha −1 a −1 in the unfertilised plots and 1.59 t ha −1 a −1 in the fertilised plots were estimated, according to Equation (3), to have been added to SOC. In the unfertilised plots, nearly 90% of the total C input was Both effects were again very highly significant (p < 0.001). The effect of undersown row crops was about half as strong as the effect of the clover-grass leys. The effect of fertilisation was the same for this subset of data as for the complete dataset. Combining all three influences for the complete data set yielded: SOCc = 1.11 + 1.7 × 10 −5 × F + 0.004 × L + 0.002 × U, R 2 = 0.46, n = 379. All effects were very highly significant (p < 0.001). R 2 indicated that these three influences explained 46% of the variation in SOCc, while 54% can be attributed to random variation and other influences (e.g., crop type, interactions).
Mean total C input, as calculated using Equation (1), was 1.98 t ha −1 a −1 in the unfertilised plots and 3.17 t ha −1 a −1 in the fertilised plots. Of this, 0.89 t ha −1 a −1 in the unfertilised plots and 1.59 t ha −1 a −1 in the fertilised plots were estimated, according to Equation (3), to have been added to SOC. In the unfertilised plots, nearly 90% of the total C input was from crop residues and roots. In the fertilised plots, 60% of the total C input came from crop residues and roots, and 33% from digestate fertilisation.
C inputs, both total and input to SOC, were highest when clover-grass was grown, due to crop residues and roots. Consequently, the calculated C input was highest for CR10, with 75% clover-grass in the rotation. Row crops such as maize and sunflower were calculated as having higher C inputs than wheat or triticale. CR8 therefore had the third-highest (after CRs 9 and 10) C inputs of the fertilised plots, particularly as the sunflower crop was undersown with clover-grass. CR3, with two years of wheat in the rotation, had the lowest C input of the fertilised plots. C inputs from digestate were highest in CRs 9 and 1, and lowest in CRs 3 and 4. The relationship between the estimated C input to SOC per crop rotation and the change in SOC for that crop rotation was very highly significant ( Figure 6) and close to the 1:1 line, even without model calibration. Note that, in particular, the data points for CRs 1a, 1b and 1c are close together, for both the fertilised and unfertilised plots.
These CRs were identical but differed in their positions within the experimental site. Their initial SOCc spanned the whole range of values for all CRs. This indicates that initial SOCc had no influence on either the calculated or measured change in SOC stocks. highest (after CRs 9 and 10) C inputs of the fertilised plots, particularly as the sunflower crop was undersown with clover-grass. CR3, with two years of wheat in the rotation, had the lowest C input of the fertilised plots. C inputs from digestate were highest in CRs 9 and 1, and lowest in CRs 3 and 4. The relationship between the estimated C input to SOC per crop rotation and the change in SOC for that crop rotation was very highly significant ( Figure 6) and close to the 1:1 line, even without model calibration. Note that, in particular, the data points for CRs 1a, 1b and 1c are close together, for both the fertilised and unfertilised plots. These CRs were identical but differed in their positions within the experimental site. Their initial SOCc spanned the whole range of values for all CRs. This indicates that initial SOCc had no influence on either the calculated or measured change in SOC stocks. Figure 6. Correlation between estimated total (sum) C input to SOC (t ha −1 ) for each CR in the period 2010-2017, and the measured change in SOC stock (t ha −1 ) of that CR between 2010 and 2017 (n = 24). Each data point is the mean of 16 plots to reduce propagation of uncertainty resulting from the subtraction of the SOC stock in 2010 from the SOC stock in 2017. Unfertilised plots are represented by grey triangles, fertilised plots by black circles. The number next to a data point is the CR. The solid line is the regression line (y = 0.61 + 0.81x, R 2 = 0.74, p < 0.001, n = 24). The dashed line is the 1:1 line. The data point marked with a square is an outlier (CR4 unfertilised plots).
Estimated C inputs to SOC varied between 5.9 and 15.5 t ha −1 depending on the CR and fertilisation, whereas the measured changes in SOC stocks ranged between 4.4 and 15.8 t ha −1 . The data indicated that the increase in SOC stocks was, on average, 1.3 t ha −1 Figure 6. Correlation between estimated total (sum) C input to SOC (t ha −1 ) for each CR in the period 2010-2017, and the measured change in SOC stock (t ha −1 ) of that CR between 2010 and 2017 (n = 24). Each data point is the mean of 16 plots to reduce propagation of uncertainty resulting from the subtraction of the SOC stock in 2010 from the SOC stock in 2017. Unfertilised plots are represented by grey triangles, fertilised plots by black circles. The number next to a data point is the CR. The solid line is the regression line (y = 0.61 + 0.81x, R 2 = 0.74, p < 0.001, n = 24). The dashed line is the 1:1 line. The data point marked with a square is an outlier (CR4 unfertilised plots).
Estimated C inputs to SOC varied between 5.9 and 15.5 t ha −1 depending on the CR and fertilisation, whereas the measured changes in SOC stocks ranged between 4.4 and 15.8 t ha −1 . The data indicated that the increase in SOC stocks was, on average, 1.3 t ha −1 less than the calculated amount of C added to the SOC stocks. The data points for both the fertilised and unfertilised plots for CR10, which had the highest proportion of clover-grass leys (three out of four years), were close to the 1:1 line, indicating that the C inputs for clover-grass were predicted correctly by the model. This was corroborated by CR9, with the second-highest proportion of clover-grass leys in the rotation (two out of four years). The data points for CR9 arranged only slightly below the 1:1 line. CRs 4 and 7, and to a lesser extent CR3, deviated most. These three CRs all had a legume/non-legume cover crop in common. However, CR8 also had a legume/non-legume cover crop but the results for this CR corresponded with the modelled values.

Discussion
The majority of studies to date investigating the effects of farm biogas production on soils have either been pot experiments or short-term conventional field experiments. Our study is based on data from a long-term field trial run according to organic principles and, as such, biogas production on the basis of clover-grass leys. It is also one of the few biogas field trials where the amount of digestate applied is based on nutrient cycling, i.e., on the estimated yield of the respective crop rotation, and which therefore represents more realistic conditions [26].
Even though the experiment had already been running for five years before the data we used in this study were collected, SOCc still increased in 96% of plots, irrespective of the CR or if plots were fertilised or not. This effect was particularly remarkable for the unfertilised plots, given that nearly all aboveground biomass had been removed for 13 years and there was no nutrient return in the form of fertilisation.
Our regression model showed that this was due to the effects of the clover-grass leys and the undersown row crops. If there were no fertilisation, clover-grass leys or undersown row crops (i.e., if F, L and U were equal to zero in Equation (6)), then our regression model would result in an SOCc of 1.11%, which is very close to the mean SOCc of 1.10% in 2010.
This effect was found across CRs and site position, and therefore any variation in soil properties. Our modelling of the C input to SOC corroborated these results and showed that an increase in SOCc based on these effects is feasible. Our results also corresponded with results from an earlier period (2010-2014) [42].
Increasing the proportion of row crops undersown with clover or clover-grass by 10% would increase SOCc by 0.02%, according to our regression model. Although only three out of 12 CRs included undersown row crops (up until 2012 CR4 also included sunflower undersown with clover-grass), the plots in these CRs still represented a wide range of differing proportions of undersown row crops. Aside from the positive effects of clover and clover-grass on SOC, undersowing row crops means more vegetation cover, which leads to more C input from roots and, in our field trial, from shoot biomass when the undersown crops were mulched. Yields of the undersown clover and clover-grass were similar in both the fertilised and unfertilised plots (data unpublished); however, yields of the main crops were lower in the unfertilised plots. Maize yields, for example, were 70% higher in the fertilised plots than in the unfertilised plots. When N is limited, crops tend to allocate more resources to the roots than aboveground biomass [43] and the root:shoot ratio increases. Nonetheless, it is likely that due to the large differences in aboveground biomass, root systems of the main crops were smaller in the unfertilised plots. The undersown crops in the unfertilised plots could therefore have produced more belowground biomass, as there was less competition for space and nutrients. Cong et al. [44], for example, found that although undersown legume-grass leys produced less aboveground biomass in unfertilised plots, there was no difference in root biomass between fertilised and unfertilised plots.
Fertilisation with digestate had the biggest effect on SOCc. In addition to the direct effect of the digestate itself on SOC, fertilisation increased yields in our trial, meaning C input from root and harvest residues was modelled as being higher in the fertilised plots. The effect of digestate fertilisation on root growth is, however, at present unclear, particularly as most studies to date have been pot experiments and data from field trials is lacking. Digestate fertilisation increased the root biomass of wheat in a field trial [45] and grass in a pot experiment [46], but in other pot experiments digestate fertilisation either had no effect [47] or inhibited root growth of grasses [48]. If both the solid and liquid phases of the digestate had been used instead of just the liquid phase as in our field trial, the effect on SOCc would have been even greater according to Equation (6), given the larger amount of C applied to the soil [49,50]. Further research would, however, be needed to ascertain if Equation (6) also applies to the solid digestate phase.
The highest inputs to SOC from crops in our field trial were when clover-grass was grown. CR10, with 75% clover-grass in the rotation, had the highest SOCc and the largest increase in SOCc between 2009 and 2017, which agrees with the finding that increasing the proportion of leys in organic CRs (also in combination with organic fertilisers) increased SOCc [28]. This increase seemed to be constant, even over a 13-year period, with an increase of 0.004% in SOCc when the proportion of clover-grass leys in the CR was increased by 1%. Given that in most rotations there was only one year of clover-grass, whereas with 75% clover-grass ley in CR10 there was no tillage in these plots for three of the four years of the CR, this indicated that tillage did not have an important influence on this effect, although it is often assumed that tillage destroys soil aggregates and exposes the SOC protected in these aggregates to mineralisation, increasing decomposition rates [51]. An important implication of the linear increase of SOCc with the proportion of clover-grass in the CR is that the length of the period between clover-grass years seemed to be of minor importance and did not cause any substantial loss of SOC. The longest break between clover-grass years was only three years in our field trial, and hence this inference is restricted to a rather short period, but given that turnover is always highest in the first few years, a low degradation rate can be expected to last longer. Freibauer et al. [52] estimated that increasing the duration of leys would lead to an increase in the potential SOC accumulation rate of 0.1-0.5 t ha −1 a −1 . Our results are consistent with this estimate. Using the average bulk density to calculate C stocks, the plots in CR10 accumulated 1.58 t ha −1 a −1 . In comparison, the plots with only one year of clover-grass ley accumulated 1.26 t ha −1 a −1 , a difference of 0.32 t ha −1 a −1 . It is important to note, however, that in CR10 ley composition changed over time, with the proportion of grass increasing, and total yield decreasing, by the third year (data unpublished). Nevertheless, we have not had any disease or pest problems in this or any of the other rotations, despite the short break between leys.
As our regressions suggest that any addition of ley to a rotation will increase SOCc, for the CRs with only 25% clover-grass leys in the rotation and no undersown row crops, it appears that the C input from the root and harvest residues from these leys was enough to increase SOCc, even without additional fertilisation. Models have estimated that C inputs of 1.5 t are enough to maintain SOC [53,54]. This is lower than the mean C input of 1.8 t ha −1 a −1 in these unfertilised plots, which could explain why SOCc did not decrease in these plots.
The correlation between the measured change in SOC and the estimated C input to SOC indicated an average respiration rate of 0.3% per year. This is lower than the mineralisation rates normally seen for arable soils [55]. Our model could therefore have underestimated the C input to SOC. The crop parameters we used were developed for conventionally-grown crops; however, it has been shown that root:shoot ratios are higher for organically-grown crops, presumably due to lower nutrient (particularly N) availability [56]. This would also explain why the difference between the measured change in SOC and the estimated C input for the unfertilised plots was smaller than for the fertilised plots, and why the measured change was higher than the C inputs in nearly half of the CRs for the unfertilised plots. However, increasing the yield-dependent parameter F RHR to take into account higher root biomass at lower yields in organic crops would only have a small effect on the estimated C input to SOC, as the parameter K RHR generally makes up a much larger proportion of the calculated C input to SOC. As K RHR is a constant independent of yield, this already allocates a larger proportion of biomass to belowground biomass at low yields.
Alternatively, it could be that the mineralisation rates in these plots were so low that barely any C was lost due to mineralisation, and therefore the C input to SOC was roughly equal to the increase in SOC. In CR10 with 75% clover-grass in the rotation, for example, in both the fertilised and unfertilised plots the measured change in SOC was greater than the estimated C input. Acharya et al. [57] found that mineralisation rates were lower in the second and third years of leys than in the first year, and that root biomass increased with ley age. This was not taken into account in our model (we used the same crop parameters for all leys), and therefore our estimated C input to SOC could have been too low for CR10.
The unfertilised plots in CR4 were an outlier in this correlation, with a much lower increase in SOC than predicted by the level of C inputs. This CR includes soybean, which is a relatively new crop in Germany. We assumed the same crop parameters for soybean as for field bean, which is a well-established crop in Germany. However, despite higher yields for soybean in our field trial than for field bean, it could be that, under German conditions, the harvest and root residues are lower. Nonetheless, reducing K RHR and the synthesis coefficient for soybean in our model, or even removing this CR from our model entirely, did not change the slope of the regression and only increased R 2 to 0.80.
Differences in yields and SOC between CRs 1a, b and c, despite identical treatment, indicate that there could be soil heterogeneity across the site. Further analysis of pH values indicated the presence of carbonates in a distinct pattern across the site. Analysis of historical records showed the presence of an old road running through the field trial. Although the 2017 soil samples were reanalysed depending on pH, this was not the case for the 2009 samples. Removing carbonates from the 2009 samples would have lowered the SOC values for these plots, and would make the difference between SOCc in 2010 and 2017 even larger. A preliminary analysis indicated, however, that it would not change the results of modelling [58]. This needs to be taken into account in the future when samples are analysed. However, this possible heterogeneity also means that our results would be valid for a wide range of initial soil conditions.
The effects of the three predictors in our regression were additive, also over time. This means our results can be applied to a broad range of situations. For example, if there are no undersown crops in the CR (i.e., U = 0 in Equation (6)), the effects of the leys on SOCc will still remain the same. Our regressions show farmers how increasing the proportion of clover-grass and/or undersown crops in their CRs will impact SOCc. The regressions can therefore also be used by governments when planning subsidy schemes for farmers with the aim of improving SOCc (for example, for C sequestration and climate change mitigation). Our additive model also means that separate subsidies could be paid for each of our three variables, so even if farmers only include one of our variables in their farming, they would still increase SOCc and could also receive a subsidy. Our results hold true for a period of at least ten years, which is a longer period than most farmers plan for and is longer than the duration of most farm subsidy schemes.
In the CANDY model, calculating the C input to SOC is also additive, where the C input from various sources is simply added together. CANDY is based on data from long-term field experiments. Despite possible sources of error, such as assuming that bulk density was the same for all plots and very simply adapting the crop parameters, our model correlated well with our measured data. This corroborates that a simple additive model based on long-term field trial data can be used to predict changes in SOCc.

Conclusions
Our long-term experiment showed that clover-grass leys have a remarkably positive effect on SOCc, increasing it by 0.004% for every year a ley is grown, even if the aboveground biomass is removed for renewable energy production and not returned as organic fertiliser. Returning the nutrients and part of the organic matter as biogas digestate increased biomass production and, in turn, SOCc even further, by 0.017% for every 1 t ha −1 of digestate C applied. When row crops were undersown with clover or clover-grass this also had a positive effect on SOCc, increasing it by 0.002% every year undersown crops were grown. Digestate fertilisation gives organic farmers the opportunity to control the timing and quantity of N fertilisation, improving yields of other crops and hence increasing C inputs to soil. Our results therefore support scenarios for sustainable farming involving the anaerobic digestion of agricultural feedstocks such as ley and crop biomass.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.