Continuous Co-Digestion of Agro-Industrial Mixtures in Laboratory Scale Expanded Granular Sludge Bed Reactors

: Anaerobic co-digestion often improves the yields and stability of single anaerobic digestion. However, ﬁnding the right substrate proportions within mixtures and corresponding optimal operating conditions using a particular reactor technology often presents a challenge. This research investigated the anaerobic digestion of three mixtures from the liquid fractions of piglet manure (PM), cow manure (CWM), starch wastewater (SWW), and sugar beet (SBT) using three 30 L expanded granular sludge bed (EGSB) reactors. The synergistic effects of two three-substrate mixtures (i.e., PM+CWM+SWW and PM+CWM+SBT) were studied using the PM+CWM mixture as a benchmark. These were used to detect the predicted synergistic interactions found in previous batch tests. The methane productivity of both three-substrate mixtures (~1.20 L CH4 /L react /d) was 2 × the productivity of the benchmark mixture (0.64 L CH4 /L react /d). Furthermore, strong indications of the predicted synergistic effects were found in the three-substrate mixtures, which were also stable due to their appropriate carbon-to-nitrogen ratio values. Moreover, the lowest averaged solid to hydraulic retention times ratio calculated for samples obtained from the top of the reactors was > 1. This conﬁrmed the superior biomass retention capacity of the studied EGSB reactors over typical reactors that have been used in agricultural biogas plants with a continuous stirred tank reactor.


Introduction
Anaerobic digestion (AD) is an efficient and suitable method for the sustainable management of bio-wastes as well as the production of biofuel [1]. It is a biological degradation process whereby biomass is converted into a mixture of gases called biogas, which consist mainly of methane and carbon dioxide, by the action of a microorganism consortium in the absence of oxygen. It is typically divided into four main stages, including hydrolysis, acidogenesis, acetogenesis, and methanogenesis [2].
However, AD is a very complex and sensitive process, which involves diverse microorganism groups that require different environmental and operational conditions [1]. For example, biomass substrate digestibility and biogas production are significantly affected by the substrate composition and chemistry, such as the carbon-to-nitrogen ratio (C/N), mineral and volatile fatty acid composition, and pH [3]. These are also affected by the operational conditions, including the hydraulic retention time (HRT), substrate loading rate, reactor temperature, and so on. Piglet manure (PM); cow manure (CWM); sugar beet (SBT); starch wastewater (SWW).

Bioreactor Setup and Operation
Three EGSB reactors with a height-diameter ratio of 3 units were employed in a continuous operation mode. The reactors were inoculated with 20 L of mesophilic inoculum with a spherical shape and dark green color. The EGSB reactors were operated at six different HRTs for 15, 10, 7, 5, 3, and 1 day(d). The HRTs were automatically altered by changing the feeding time of the pump at a constant flow. The HRT was calculated by Equation (1) [19,25].
where the HRT, volume of the reactor (V R ), and influent volumetric flow are in day(s), m 3 , and m 3 /day units, respectively. The recirculation pump was continuously working at an up-flow velocity of 5 m/h. Each reactor was connected to a 100 L tank that was kept under a nitrogen atmosphere and temperature of 4 • C to prevent premature aerobic degradation. The scheme for a single reactor is shown in Figure 1.
All three reactors were operated under mesophilic conditions with temperatures between 37 and 40 • C and pH values close to 8 by regulating the feed of the reactor, which is mainly possible due to the buffer capacity of the manures [26,27]. The procedure for the settings and monitoring of the continuous operation was as described in reference [28]. The measured and set variables are summarized in Table 2. Other relevant values were calculated from registered variables such as the organic loading rate (OLR) (kg COD /m 3 /d), methane productivity (MPR) (L CH4 /L reactor /d)), methane yield (MY) (L CH4 /kg VS ), removal efficiencies of chemical oxygen demand ( where the HRT, volume of the reactor (VR), and influent volumetric flow are in day(s), m 3 , and m 3 /day units, respectively. The recirculation pump was continuously working at an up-flow velocity of 5 m/h. Each reactor was connected to a 100 L tank that was kept under a nitrogen atmosphere and temperature of 4 °C to prevent premature aerobic degradation. The scheme for a single reactor is shown in Figure 1.
All three reactors were operated under mesophilic conditions with temperatures between 37 and 40 °C and pH values close to 8 by regulating the feed of the reactor, which is mainly possible due to the buffer capacity of the manures [26,27]. The procedure for the settings and monitoring of the continuous operation was as described in reference [28]. The measured and set variables are summarized in Table 2. Other relevant values were calculated from registered variables such as the organic loading rate (OLR) (kgCOD/m 3 /d), methane productivity (MPR) (LCH4/Lreactor/d)), methane yield (MY) (LCH4/kgVS), removal efficiencies of chemical oxygen demand (ŋCOD) (%), and biological oxygen demand on the fifth day (ŋBOD5) (%).
x Chemical oxygen demand (mg O2 /L) x x Biochemical oxygen demand at 5th day (mg O2 /L) x x Loading rate per unit volume (kg COD /m 3 ·d) x Hydraulic residence time (d) x pH value (-) x Gas composition in volume fractions (%, ppm) x Ratio of volatile organic acids to total inorganic carbon (-) x x: measured in the coresponded position.

Data Cleaning and Analysis
Data cleaning was performed with the aim of obtaining a data set which did not contain obvious failures, start-up periods, or clear mistakes in an operation. The data cleaning was performed using the three main criteria as described below.
The chemical oxygen demand removal is ≥0.

Overview of Each Reactor's Operation
A comparison of the different operation points of a given reactor was made using the variables of MY, MPR, Appl. Sci. 2022, 12,  ). An analysis of the practical operation of each reactor was completed using this information and the complementary information on the mixtures involved. Box and scatter plots were employed to visualize each reactor's operation.

Principal Component Analysis
PCA is an adaptive exploratory method which can be used on numerical data of various types. From a mathematical point of view, principal components are linear combinations of original variables, making them orthogonal to each other [29,30]. This method increases the interpretability of the data and at the same time minimizes information loss. For each reactor, a new data set was created using the average values of all the operation points for the above for all four response variables. The new datasets were used in a principal component analysis (PCA) to compare the reactors in terms of operation points as part of a multivariate analysis. Up to five components were acceptable and three components were desirable. The goal of the PCA was to rotate the data into an axis system where the greatest amount of variance was captured in a small number of dimensions [31].
The PCA involved the calculation of the eigenvectors and eigenvalues of a sample covariance or correlation matrix. Furthermore, the calculation of the principal components was carried out using a singular value decomposition (SVD) [32]. PCA was also employed for outlier detection due to its robustness [32,33].

SRT/HRT
A high-rate reactor such as an ESGB reactor can decouple HRT and SRT, thereby increasing the residence time of a biomass element within the reactor [15,34,35]. One of the main selection criteria for a reactor is a high SRT/HRT ratio, which prevents the washout of slow-growing methanogens [15]. The sludge age (SRT) in d is given by Equation (2).

SRT =
Mass of sludge in reactor Mass of sludge wasted per day (d) (2) If a steady-state condition was assumed, Equation (2) can be written as Equation (3).
where x i , V R , Q eff , and X eff are the viable biomass concentration inside the reactor (kg VS / kg FM ), volume of the bioreactor (m 3 ), effluent flow out (L/d) of the reactor, and viable biomass concentration in the effluent (kg VS /kg FM ), respectively. Since the input and output flows were equal (steady-state condition), Equation (3) was transformed to (4).
The ratio SRT/HRT values were calculated for all three reactors using a biomass, which was sampled from the top of each reactor where biomass concentration was lowest. The ideal SRT/HRT ratio should be >3 [15,17,19].

Characterization of Synergistic Effects
The synergistic effects of the three-substrate mixtures (PM+CWM+SWW and PM+ CWM+SBT) were compared using the two-substrate mixture (PM+CWM) as a benchmark. Since it was not possible to operate a single digestion of each substrate, the hypothesis employed used Equation (5) (5) where MY MAX_TM and MY MAX_DM are the maximum MY of the three-and two-substrate mixtures, respectively. Since the calculated ratios were based on the yield of the two-substrate benchmark mixture, the ratio for the PM+CWM mixture was equal to 1. As for the batch data results from Regalado et al. [24] and some complementary unpublished data, the MY MAX value corresponding to a maximum value of each operating point was used. In addition, a comparison between the MY MAX in the continuous and batch operation processes was made for each mixture.

Characterization of Hydraulic Behaviors
The hydrodynamics of the anaerobic reactor was studied because they significantly influence the rates of biological reactions. They particularly affect the rates of mass transfer and the distribution of reactions along a reactor, both of which determine a reactor's overall performance [20,36]. The amount of mixing in a reactor also determines the performance of a reactor; therefore, to describe the real behavior of a reactor, the influence of mixing on the mass balance equation must be specified correctly [37]. In this study, the hydrodynamics were characterized by the non-dimensional numbers given by Peclet and Reynolds.
The mixing intensity of the fluid within a reactor is well described by the axial Peclet number (Pe axial ) (see Equation (6)).
where V up , H, and D A are the up-flow velocity (m/h), bioreactor height (m), and axial dispersion coefficient (m 2 /h), respectively. When D A → ∞, the value of Pe axial becomes 0 since Pe axial is an inverse function of D A . Consequently, the system will operate as a plug-flow reactor since there is no mixing in the axial direction. On the other hand, when D A → 0, the system will behave as a complete mixture reactor [19]. Various transfer functions have been proposed to estimate the dispersion from either the Reynolds number or a flow velocity [38][39][40][41][42]. Here, we used an approach described by Equations (7) and (8), which assessed D A as a function of flow distance D A = 1.03·V up 1.11 ·0.009 n j (7) where n j , z, and H are the values of the normalized height, axial position (m), and height (m) of the bioreactor, respectively. The amount of turbulence is characterized by the Reynolds number (Re) and is given by Equation (9) [19,43]. Where V up , d, µ w , ρ w , and υ w are the up-flow velocity (m/h), bioreactor diameter (m), dynamic viscosity (Pa·s), density (kg/m 3 ), and kinematic viscosity (m 2 /s), respectively, Reynolds describes a relationship between inertial to viscous forces [42]. Equation 9 is the most widely used; however, it exits variations of the Reynolds number around noncircular conduits, packed beds, and mixing impellers.
Turbulence, meanwhile, is the rotational and three-dimensional chaotic movement in all directions of flowing elements, where the resulting net flow is unidirectional. The rapid mixing associated with turbulence enhances the momentum, heat, and mass transfer processes. The intervals of Reynolds include Re < 2300, 2300 < Re < 4000, and Re > 4000, which correspond to laminar, transient, and turbulent regimen, respectively. However, a typical turbulent regimen truly manifests itself from values of Re > 10,000.

Stover-Kincannon Model
The MPRs of the reactors were modeled using a variation of the Stover-Kincannon model for an anaerobic filter reactor (Equation (10)), which was proposed and implemented by Yu et al., Verma et al., Jafarzadeh et al. [44][45][46].
where MPR, MPR max , and M B are the methane productivity (L CH4 /L react /d), maximum MPR (L CH4 /L react /d), and constant (kg COD /m 3 ·d), respectively. OLR is the organic loading rate (kg COD /m 3 ·d). A non-linear regression procedure was employed using the calculated clean averaged data of all reactors. To identify similarities and differences in the kinetic behavior of all possible combinations, the averaged data of reactors 1, 2, and 3 were arranged to have a total of seven datasets. The goodness of fit was measured by a root-meansquare error (RMSE). For the most meaningful dataset(s), a simple regression analysis was performed for MY and Once the significant models were identified, their dependencies were plotted with OLR to perform a graphical optimization. In a graph, the ordinates represented the values of the individual variable divided by their maximum measured value (V i /V Max ), which was expressed in%. Thus, the ordinates represented values between 0 and 100% for each plotted variable.

Reactors' Operation Overview
Data cleaning was performed according to the set criteria in Section 2.3. The twosubstrate mixture of PM+CWM (pellets 1) had a BMP ∞ of 342.83 L CH4 /kg VS . Meanwhile, the BMP ∞ values for the three-substrate mixtures of pellets 2 and 3 were 534.21 and 530.28 L CH4 /kg VS , respectively [24,47]. After data cleaning, the resulting data sets have sizes of 162, 181, and 203 instances for reactors 1, 2, and 3 respectively. To accurately characterize the performance of a reactor, four employed output variables, including MPR (L CH4 /L reactor /d), MY (L CH4 /kg VS ),

HRT
Hydraulic retention time AcoD Anaerobic co-digestion SRT  [28,48]. The overview of reactor 1 is shown in Figure 2. The MPR suffered a sudden drop at the operating point of 5 d. This interrupted the upward trend that was observed from HRT of 15 to 7 d. From 3 to 1 d, a sustained increase in the MPR was observed, with a value that was almost 2× the second-highest average value observed at 7 d. This suggested that a punctual failure occurred at 5 d, which was not due to reaching an operational HRT limit. MY had the highest mean value of 272 L CH4 /kg VS at 15 d. After 15 d, the values significantly decreased and a similar sudden drop in MPR was observed at a HRT of 5 d. However, MY did not experience a significant recovery after the inhibition, unlike MPR. Thus, by considering both variables of MPR and MY simultaneously, the operation can be divided into two main stages: before and after inhibition. There is a noteworthy difference between these two stages. The former reached considerably higher yields than the latter; however, similar values of MPR were found in both stages.
The removal efficiencies in these two stages were not as evident. BOD 5 removal efficiency values noticeably dropped at 10 and 1 d. The low BOD 5 removal values may explain the drop in MY at 10 and 1 d in the previous operation point. This was most likely due to low reaction completions [49]. A minimum average value for the chemical oxygen demand (COD) removal efficiency was observed at 5 d. However, the trend followed by the average

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

Conflicts of Interest:
The authors declare no conflict of interest.

Reactors' Operation Overview
Data cleaning was performed according to the set criteria in Section 2.3. The twosubstrate mixture of PM+CWM (pellets 1) had a BMP∞ of 342.83 LCH4/kgVS. Meanwhile, the BMP∞ values for the three-substrate mixtures of pellets 2 and 3 were 534.21 and 530.28 LCH4/kgVS, respectively [24,47]. After data cleaning, the resulting data sets have sizes of 162, 181, and 203 instances for reactors 1, 2, and 3 respectively. To accurately characterize the performance of a reactor, four employed output variables, including MPR (LCH4/Lreactor/d), MY (LCH4/kgVS), ŋCOD, and ŋBOD5 (%), were used [28,48]. The overview of reactor 1 is shown in Figure 2. The MPR suffered a sudden drop at the operating point of 5 d. This interrupted the upward trend that was observed from HRT of 15 to 7 d. From 3 to 1 d, a sustained increase in the MPR was observed, with a value that was almost 2x the secondhighest average value observed at 7 d. This suggested that a punctual failure occurred at 5 d, which was not due to reaching an operational HRT limit. MY had the highest mean value of 272 LCH4/kgVS at 15 d. After 15 d, the values significantly decreased and a similar sudden drop in MPR was observed at a HRT of 5 d. However, MY did not experience a significant recovery after the inhibition, unlike MPR. Thus, by considering both variables of MPR and MY simultaneously, the operation can be divided into two main stages: before and after inhibition. There is a noteworthy difference between these two stages. The former reached considerably higher yields than the latter; however, similar values of MPR were found in both stages.
The removal efficiencies in these two stages were not as evident. BOD5 removal efficiency values noticeably dropped at 10 and 1 d. The low BOD5 removal values may explain the drop in MY at 10 and 1 d in the previous operation point. This was most likely due to low reaction completions [49]. A minimum average value for the chemical oxygen demand (COD) removal efficiency was observed at 5 d. However, the trend followed by the average ŋCOD values had a smaller variation compared to both MPR and MY. Furthermore, both BOD5 and removal efficiencies behaved differently since the calculated BOD5/COD ratios were fluctuating.

Reactors' Operation Overview
Data cleaning was performed according to the set criteria in Section 2.3. The twosubstrate mixture of PM+CWM (pellets 1) had a BMP∞ of 342.83 LCH4/kgVS. Meanwhile, the BMP∞ values for the three-substrate mixtures of pellets 2 and 3 were 534.21 and 530.28 LCH4/kgVS, respectively [24,47]. After data cleaning, the resulting data sets have sizes of 162, 181, and 203 instances for reactors 1, 2, and 3 respectively. To accurately characterize the performance of a reactor, four employed output variables, including MPR (LCH4/Lreactor/d), MY (LCH4/kgVS), ŋCOD, and ŋBOD5 (%), were used [28,48]. The overview of reactor 1 is shown in Figure 2. The MPR suffered a sudden drop at the operating point of 5 d. This interrupted the upward trend that was observed from HRT of 15 to 7 d. From 3 to 1 d, a sustained increase in the MPR was observed, with a value that was almost 2x the secondhighest average value observed at 7 d. This suggested that a punctual failure occurred at 5 d, which was not due to reaching an operational HRT limit. MY had the highest mean value of 272 LCH4/kgVS at 15 d. After 15 d, the values significantly decreased and a similar sudden drop in MPR was observed at a HRT of 5 d. However, MY did not experience a significant recovery after the inhibition, unlike MPR. Thus, by considering both variables of MPR and MY simultaneously, the operation can be divided into two main stages: before and after inhibition. There is a noteworthy difference between these two stages. The former reached considerably higher yields than the latter; however, similar values of MPR were found in both stages.
The removal efficiencies in these two stages were not as evident. BOD5 removal efficiency values noticeably dropped at 10 and 1 d. The low BOD5 removal values may explain the drop in MY at 10 and 1 d in the previous operation point. This was most likely due to low reaction completions [49]. A minimum average value for the chemical oxygen demand (COD) removal efficiency was observed at 5 d. However, the trend followed by the average ŋCOD values had a smaller variation compared to both MPR and MY. Furthermore, both BOD5 and removal efficiencies behaved differently since the calculated BOD5/COD ratios were fluctuating. To identify the potential causes of inhibition, the control variables of OLR and VOA/TIC were employed. The results are summarized in Figure 3. The high OLR values observed were not an indicative cause of inhibition. In particular, the OLR value was not very high at 5 d. Moreover, the system did not run at very high OLR values despite a constant decrease in the HRT values. Therefore, the COD values in the inflow suffered sizable fluctuations, as shown in Figure 4.
OLR fluctuations were not strongly correlated with the VOA/TIC, which implies that a failure was not due to a system overload caused by overfeeding. Instead, they seemed to be more connected to the quality of the feed (Figure 4). However, the VOA/TIC results showed that an acid accumulation occurred at 10 d. In these results, the values decreased sharply and approached zero at the operating point of 5 d. Hence, two possibilities were weighted, such as the non-failure of a system due to OAs accumulation and a system delayed response due to VOA accumulation at 10 d. In the former, there was no sign of strong inhibition before 5 d and a lack of VOAs in the system. This were supported by the close to zero VOA/TIC values. In the latter, the system had a delayed response to VOA accumulation at 10 d, which seemed more unlikely, given how large the delay had to be. In addition, the MPR values increased from 10 to 7 d, while the MY values were practically the same. Therefore, either the acetogenesis was the limiting-rate step or the quality of the feed was very low. To identify the potential causes of inhibition, the control variables of OLR and VOA/TIC were employed. The results are summarized in Figure 3. The high OLR values observed were not an indicative cause of inhibition. In particular, the OLR value was not very high at 5 d. Moreover, the system did not run at very high OLR values despite a constant decrease in the HRT values. Therefore, the COD values in the inflow suffered sizable fluctuations, as shown in Figure 4.
OLR fluctuations were not strongly correlated with the VOA/TIC, which implies that a failure was not due to a system overload caused by overfeeding. Instead, they seemed to be more connected to the quality of the feed (Figure 4). However, the VOA/TIC results showed that an acid accumulation occurred at 10 d. In these results, the values decreased sharply and approached zero at the operating point of 5 d. Hence, two possibilities were weighted, such as the non-failure of a system due to OAs accumulation and a system delayed response due to VOA accumulation at 10 d. In the former, there was no sign of strong inhibition before 5 d and a lack of VOAs in the system. This were supported by the close to zero VOA/TIC values. In the latter, the system had a delayed response to VOA accumulation at 10 d, which seemed more unlikely, given how large the delay had to be. In addition, the MPR values increased from 10 to 7 d, while the MY values were practically the same. Therefore, either the acetogenesis was the limiting-rate step or the quality of the feed was very low.  COD in and VOA s_in values were determined for the feed to investigate whether the quality of the feed was responsible for the inhibition. The results are summarized in Figure 4. All the targeted acids were found except for valeric and caproic acids. A rapid decline in the total acetic acid concentration and equivalent occurred after 10 d. Likewise, the concentration of the acids was almost zero at 5 d. This behavior was consistent with the inhibition observed at the operating point of 5 d. Thus, the occurrence of a failure due to the lack of VOAs in the feed, which caused very low VOA/TIC values, was accepted.
The observed COD in fluctuations explained the behavior of the OLR with the gradual reduction of HRTs. It was expected that intensive variables, such as COD in , VOA concentrations, and the BOD 5 /COD ratio, would show stable behavior. However, these variables fluctuated due to the lack of proper mixing in the feeding tank. Continuous mixing was not done during operations, although the mixtures were vigorously mixed in the tank during preparation. This induced the settling of particulates and instability of the feed, which caused the first excess of VOAs observed at 10 d. The concentration of VOAs was approximately zero at 5 d, which suggested the existence of a substrate limitation on the system. This limitation substantially influenced the operation since the tank had to be refilled several times, causing variations in the preparation. Therefore, heterogeneities in the composition of substrate mixtures during a year of operation are expected due to seasonal behavior [50]. Appl. Sci. 2022, 12, x FOR PEER REVIEW 10 of 23 CODin and VOAs_in values were determined for the feed to investigate whether the quality of the feed was responsible for the inhibition. The results are summarized in Figure 4. All the targeted acids were found except for valeric and caproic acids. A rapid decline in the total acetic acid concentration and equivalent occurred after 10 d. Likewise, the concentration of the acids was almost zero at 5 d. This behavior was consistent with the inhibition observed at the operating point of 5 d. Thus, the occurrence of a failure due to the lack of VOAs in the feed, which caused very low VOA/TIC values, was accepted.
The observed CODin fluctuations explained the behavior of the OLR with the gradual reduction of HRTs. It was expected that intensive variables, such as CODin, VOA concentrations, and the BOD5/COD ratio, would show stable behavior. However, these variables fluctuated due to the lack of proper mixing in the feeding tank. Continuous mixing was not done during operations, although the mixtures were vigorously mixed in the tank during preparation. This induced the settling of particulates and instability of the feed, which caused the first excess of VOAs observed at 10 d. The concentration of VOAs was approximately zero at 5 d, which suggested the existence of a substrate limitation on the system. This limitation substantially influenced the operation since the tank had to be refilled several times, causing variations in the preparation. Therefore, heterogeneities in the composition of substrate mixtures during a year of operation are expected due to seasonal behavior [50]. Similar trends to reactor 1 were observed for reactors 2 and 3. For reactor 2 in particular, the inhibition was the least abrupt due to a lower dry matter (DM) content (i.e., almost 80% of SWW on a fresh matter basis), which reduced the effect of the seasonal behavior observed in the manures (Table 1). A summary of the averaged behavior for a specific operation point for each reactor is found in Table 3, which uses a three-color scale by column. The colors were ordered red, yellow, and green to show the increase from lower to higher values. The intensity of each color was determined by its proximity to the lowest, middle, or highest value. Similar trends to reactor 1 were observed for reactors 2 and 3. For reactor 2 in particular, the inhibition was the least abrupt due to a lower dry matter (DM) content (i.e., almost 80% of SWW on a fresh matter basis), which reduced the effect of the seasonal behavior observed in the manures (Table 1). A summary of the averaged behavior for a specific operation point for each reactor is found in Table 3, which uses a three-color scale by column. The colors were ordered red, yellow, and green to show the increase from lower to higher values. The intensity of each color was determined by its proximity to the lowest, middle, or highest value.  The colors are ordered red, yellow, and green from lower to higher values. R is the reactor.
The colors are ordered red, yellow, and green from lower to higher values. R is the reactor. Computer flow dynamics COD removal efficiency observed in reactor 1 was much smaller than in reactors 2 and 3. The same trend was observed with MY. The COD removal efficiency has been interpreted as a degree of reaction completion [19,41,49]. Thus, a strong relationship between COD removal and reaction completion was expected. Meanwhile, a similar relationship should have been found between MY and ; however, due to fluctuations in the BOD 5 /COD, no clear visual correlation could be established. Furthermore, the drop in the MYs of all three reactors at HRT 5 d was linked with the lack of organic dry matter (ODM) or VOA in the feed, as previously shown in Figure 4.
The MYs in reactors 2 and 3 were also significantly higher than those in reactor 1. The maximum MY of reactors 2 and 3 were similar at 7 and 10 d, respectively. The improvement observed thanks to the addition of a third substrate to the PM+CWM mixture was probably related to a higher C/N ratio. The C/N ratio balance in feedstocks was significant for the stable operation of AD. Substrates with high C/N ratios have a poor buffering capacity; therefore, nitrogen will be consumed rapidly by methanogens to meet their protein requirements. This results in low methane production and produces excess VOAs during fermentation. In typical feedstocks with a low C/N ratio, nitrogen has been found to accumulate in the form of ammonia, which inhibits the methanogens and prevents methane production [4,51].
The lowest recommended limit for C/N is 20 [11]; thus, a value of 15 was sufficient for our purpose [13]. The values of C/N in Table 1 are between 15 and 20 for the threesubstrate mixtures. Meanwhile, the value was <15 for the two-substrate mixture. This supported our finding that the three-substrate mixtures were more stable and produced more methane. Nonetheless, Lissens et al. [52] have affirmed that substrates with a C/N ratio <10 can support a stable process; however, they require a multistage system to avoid reactor overloading.
The maximum values observed for the COD and BOD 5 removal efficiencies, as well as the MY for reactors 2 and 3, were at HRTs of either 10 or 7 d. The same operation interval was observed by Cruz-Salomón et al. [19] and regarded as the optimal operation interval for EGSB reactors.
All three reactors reached their maximum MPR at 1 d. The observed stable operation of the EGSB reactors at HRTs of 3 and 1 d presented some novelty in our operation with the agricultural substrates. Castrillón Cano et al. [53] were able to operate reactors at HRTs of as low as 8 h; however, they only used a 3.4 L effective volume to perform their residence time distribution (RTD) experiments with water in the presence and absence of biomass. In another study, Dereli [54] effectively operated a full-scale EGSB of 1200 m 3 at an average HRT of 7 d for the treatment and digestion of confectionery industry wastewater. Meanwhile, Cruz-Salomón et al. [25] performed continuous tests with a 3.3 L EGSB reactor with a HRT of between 3 and 9 d for the treatment of coffee processing wastewater. In addition, Rico et al. [55] operated a UASB reactor with an external settler and effluent recycling for alkalinity supplementation for the co-digestion of cheese whey and the liquid fraction of dairy manure. Under a constant HRT of 2.2 d, their system demonstrated a stable operation with up to 75% cheese whey fraction in the feed. This was with an applied OLR of 19 kg COD m −3 d −1 , obtaining a   COD and MPR values were 95.1% and 9.5 m 3 CH4 m −3 d −1 , respectively. Therefore, we conclude that there is a novelty in our successfully operated EGSB reactors for AD from the agricultural substrates in mixtures 1 and 3 at small HRTs. Notably, Rico et al. [55] suggested that a manipulation of the mixture proportion at a constant HRT can also lead to improvements in terms of both stability and MPR.

PCA
PCA was conducted using the data shown in Table 3. The results are shown in Figure 5.
The score of Figure 5a shows the distribution of the data in the reactor number combining the shapes and colors shown in the figure legend. The HRT is shown above each point. The axes of the graph were created by a linear combination of the variables. This is represented in Figure 5c with MY and COD removal being the most influential variables in the x and y-axis. respectively. Meanwhile, Figure 5b shows the combination of Figure 5a,c. The red and blue points are the variable and data points, respectively. The x-axis was the most significant since it explained most of the variability of the data (e.g., up to 99%). PCA has been used for reducing the dimensionality of large datasets [29,30,56]. However, since the dataset employed was small, PCA was used for descriptive purposes only.

PCA
PCA was conducted using the data shown in Table 3. The results are shown in Figure  5. The score of Figure 5a shows the distribution of the data in the reactor number combining the shapes and colors shown in the figure legend. The HRT is shown above each point. The axes of the graph were created by a linear combination of the variables. This is represented in Figure 5c with MY and COD removal being the most influential variables in the x and y-axis. respectively. Meanwhile, Figure 5b shows the combination of Figure 5a,c. The red and blue points are the variable and data points, respectively. The x-axis was the most significant since it explained most of the variability of the data (e.g., up to 99%). PCA has been used for reducing the dimensionality of large datasets [29,30,56]. However, since the dataset employed was small, PCA was used for descriptive purposes only. The MPR was not able to sufficiently describe the variability of the data since the most influential variables were found between the two external ellipses in Figure 5c. MY was the most important and efficient according to the PCA results within the low-right quarter of the ellipse in Figure 5a. The green diagonal line in both Figure 5a,b represent the difference between the efficient and non-efficient operations. Hence, the best points were 7, 10, and 7 d for reactors 1, 2, and 3, respectively. These points reached their high values simultaneously in all four variables ( Table 3). The operation of reactor 1 reached a comparable efficiency to reactors 2 and 3 at a HRT of 15 d only. Nonetheless, a more efficient operation at lower HRTs was possible with the three-substrate mixtures.
Most of the operating points for HRTs at <7 d were considered inefficient, except for reactor 2 with a HRT of 3 d. Therefore, the HRTs at <7 d were generally feasible; however, they are not recommended due to their low efficiencies.

SRT/HRT
EGSBs can potentially reach much lower HRTs than CSTR reactors, which are typically used in agricultural biogas plants. In these reactors, the HRTs are equal to SRTs. However, an ESGB reactor can decouple the retention times by increasing the residence time of the biomass within the reactor [15,34,35]. The SRT/HRT values were calculated by equation 4 for each reactor using a biomass sample from the top of a reactor, where the The MPR was not able to sufficiently describe the variability of the data since the most influential variables were found between the two external ellipses in Figure 5c. MY was the most important and efficient according to the PCA results within the low-right quarter of the ellipse in Figure 5a. The green diagonal line in both Figure 5a,b represent the difference between the efficient and non-efficient operations. Hence, the best points were 7, 10, and 7 d for reactors 1, 2, and 3, respectively. These points reached their high values simultaneously in all four variables ( Table 3). The operation of reactor 1 reached a comparable efficiency to reactors 2 and 3 at a HRT of 15 d only. Nonetheless, a more efficient operation at lower HRTs was possible with the three-substrate mixtures.
Most of the operating points for HRTs at <7 d were considered inefficient, except for reactor 2 with a HRT of 3 d. Therefore, the HRTs at <7 d were generally feasible; however, they are not recommended due to their low efficiencies.

SRT/HRT
EGSBs can potentially reach much lower HRTs than CSTR reactors, which are typically used in agricultural biogas plants. In these reactors, the HRTs are equal to SRTs. However, an ESGB reactor can decouple the retention times by increasing the residence time of the biomass within the reactor [15,34,35]. The SRT/HRT values were calculated by equation 4 for each reactor using a biomass sample from the top of a reactor, where the biomass concentration was lowest. A steady-state condition was assumed in the calculations and the results are shown in Table 4. tions and the results are shown in Table 4. The ideal SRT/HRT ratio should be >3 [17]; however, this was far from being accomplished by sampling from the top of a reactor. In all cases, the averaged SRT/HRT was >1. Nevertheless, in three instances (one from each reactor), SRT/HTR ratios smaller than 1 were calculated. This demonstrates that even by sampling from where the biomass concentration is lowest, an average EGSB reactor can retain a better biomass than a typical CSTR. Minimal washout was observed in all three reactors. Unfortunately, data was not collected between the HRTs of 10 and 7 d. The red and green highlighted numbers were the lowest and highest values in the columns, respectively.

Characterization of Synergistic Effects
To study the possible synergistic effects suggested by the interactions identified as acute effects in Regalado et al. [24], the ratios between the maximum MYs in the continuous operation and batch validation tests were compared using equation 5. The results are shown in Table 5. Table 5. Methane yields and ratios based on piglet and cow manure yield in the batch and continuous tests.

Maximum Methane Yield in Continuous
Tests (LCH4/kgVS) The MY ratio of the mixture with SWW was almost the same in both scales. The relative difference in the ratios obtained for the mixture with SBT was 7.01%. This confirmed the acute effects of adding a third substrate to the two-substrate mixture. The third substrate provided the same boost in the continuous stage and batch tests.

Continuous to Batch
Also, the methane yield ratio obtained during the transfer from batch tests to the continuous stage was between 0.72 and 0.78. It was expected that the obtained MY from the continuous stage would be smaller than the ultimate biomethane potential from the batch test described by Weinrich and Nelles [11]. Similar intervals in the continuous stage to batch tests methane yields ratio have been identified in the literature. For example, Mahnert et al. [57] obtained ratios from 0.73 to 0.8 from the use of different grass species. Obiukwu and Nwafor [58] reached a ratio of 0.81 from the use of grape pomace. Meanwhile, Chowdhury and Fulford [59] used the mesophilic digestion of cattle dung in both batch and semi-continuous digestion with four reactors and six semi-continuous reactors, The red and green highlighted numbers were the lowest and highest values in the columns, respectively.
The ideal SRT/HRT ratio should be >3 [17]; however, this was far from being accomplished by sampling from the top of a reactor. In all cases, the averaged SRT/HRT was >1. Nevertheless, in three instances (one from each reactor), SRT/HTR ratios smaller than 1 were calculated. This demonstrates that even by sampling from where the biomass concentration is lowest, an average EGSB reactor can retain a better biomass than a typical CSTR. Minimal washout was observed in all three reactors. Unfortunately, data was not collected between the HRTs of 10 and 7 d.

Characterization of Synergistic Effects
To study the possible synergistic effects suggested by the interactions identified as acute effects in Regalado et al. [24], the ratios between the maximum MYs in the continuous operation and batch validation tests were compared using equation 5. The results are shown in Table 5. Table 5. Methane yields and ratios based on piglet and cow manure yield in the batch and continuous tests.

Maximum Methane Yield in Continuous
Tests (L CH4 /kg VS )

Methane Yield Predicted by the Model in Batch Tests
(L CH4 /kg VS ) The MY ratio of the mixture with SWW was almost the same in both scales. The relative difference in the ratios obtained for the mixture with SBT was 7.01%. This confirmed the acute effects of adding a third substrate to the two-substrate mixture. The third substrate provided the same boost in the continuous stage and batch tests.

Continuous to
Also, the methane yield ratio obtained during the transfer from batch tests to the continuous stage was between 0.72 and 0.78. It was expected that the obtained MY from the continuous stage would be smaller than the ultimate biomethane potential from the batch test described by Weinrich and Nelles [11]. Similar intervals in the continuous stage to batch tests methane yields ratio have been identified in the literature. For example, Mahnert et al. [57] obtained ratios from 0.73 to 0.8 from the use of different grass species. Obiukwu and Nwafor [58] reached a ratio of 0.81 from the use of grape pomace. Meanwhile, Chowdhury and Fulford [59] used the mesophilic digestion of cattle dung in both batch and semi-continuous digestion with four reactors and six semi-continuous reactors, respectively.
Their results showed higher rates in the semi-continuous operation; however, biogas yields were lower compared to the batch test. Their batch tests reached 67% COD efficiencies, which was lower than the results in this study. In addition, Holliger et al. [60] suggested that BMPs can be used to estimate biogas production at full scale; however, the BMP value should be multiplied by a factor of 0.8-0.9 to avoid overestimation.

Hydraulic Analysis
The results for both hydraulic parameters of the reactors are shown in Table 6. The Pe axial results showed values that were very close to 0, even in the separation zone. This indicated a flow pattern that was close to a completely mixed system [61]. The value of the Reynolds was also very similar to the one obtained by Brito and Melo [36]. These authors fitted an EGSB reactor to a CSTR with a characteristic coefficient of determination of 0.92. The inclusion of a short circuit increased the coefficient of determination to 0.95. Therefore, a CSTR model for simplicity was accepted and successfully used in their mass balance equations. Similarly, López and Borzacconi [62] assumed a CSTR behavior based on a high recirculation ratio and expansion of a bed. The combination of these two effects resulted in the significant mixing of the liquid and solid phases, as well as uniform gas production. However, their mass balance equation for the biomass included a washout effect, which was attributed to the high up-flow velocity at which the reactor was operated.
Nevertheless, the relative increase of the Peclet's value from one zone to another was considered significant. Consequently, due to the different behaviors of the zones in the reactor, the zones can be modeled as different reactors in series as described by Gleyce et al. [20]. Gleyce et al. [20] divided a reactor into two major zones, i.e., the separator and the reactor tubes. The reactor could be modeled either as two plug-flow reactors in series or five CSTRs with three separators and two tubes with coefficients of determination of 0.94 and 0.95, respectively.

Modeling of Reactors
The Stover-Kincannon model was used to fit the five different combinations of data from the reactors. The averaged values were employed and the size of a combination was up to 18 points. The results are summarized in Table 7. The fit for R2 was very good and the best among all datasets. The datasets for R1 and R3 had the worst fit, which suggested that the largest difference in the kinetic behavior among all possible combinations existed in these reactors. The RMSE values for R1 and R2 were also rather large, which meant that there were significant differences in the behaviors of these reactors. The goodness of fit in R1 and R3 were equal; however, R3 could produce a maximum amount of methane that was more than double the daily amount of methane from R1. R3 was also able to handle a larger OLR, which was related to the intrinsic properties of the substrate mixtures. In addition, the model predicts that R2 was far from reaching its maximum in terms of the production that can be handled by its largest OLR. The datasets for R2 and R3 were of special interest. Since the mixtures in these reactors showed synergistic effects and they seemed to behave similarly, we examined if the mixtures followed similar kinetics. We found that the difference among them was moderate (R2, R3); however, the differences between R1 and R2 (R1, R2) or R1 and R3 (R1, R3) were larger. The fits and measured data are shown in Figure 6. By following both estimated models (red and violet lines), it was noticed that significant differences existed at low OLRs, with these being much smaller at higher OLRs. The fluctuations in R3 between 4.5 and 6.5 gCOD/L were most likely the main cause of the misfit, which was observed in the green but not in the red line. The performance of R3 (green line) was closely related to mixture preparation and the degree of mixing in the feeding tank since mixture 3 had the highest DM content, contrary to the smooth behavior of R2 (light blue line). Appl. Sci. 2022, 12,   The datasets for R2 and R3 were of special interest. Since the mixtures in these reactors showed synergistic effects and they seemed to behave similarly, we examined if the mixtures followed similar kinetics. We found that the difference among them was moderate (R2, R3); however, the differences between R1 and R2 (R1, R2) or R1 and R3 (R1, R3) were larger. The fits and measured data are shown in Figure 6. By following both estimated models (red and violet lines), it was noticed that significant differences existed at low OLRs, with these being much smaller at higher OLRs. The fluctuations in R3 between 4.5 and 6.5 gCOD/L were most likely the main cause of the misfit, which was observed in the green but not in the red line. The performance of R3 (green line) was closely related to mixture preparation and the degree of mixing in the feeding tank since mixture 3 had the highest DM content, contrary to the smooth behavior of R2 (light blue line). Therefore, we decided to work on both mixtures individually, given that both models were not able to converge in most of the working intervals. Since working at low HRTs usually reduces the MY [44], we took into account the other response variables in order to establish an optimal operational OLR. Hence, empirical models of MY and ŋCOD versus OLR were also fitted. The fits for reactors 2 and 3 are shown in Figure 7a and 7b, respectively. Therefore, we decided to work on both mixtures individually, given that both models were not able to converge in most of the working intervals. Since working at low HRTs usually reduces the MY [44], we took into account the other response variables in order to establish an optimal operational OLR. Hence, empirical models of MY and  No significant fit was found for ŋCOD, as per the pre-established criterion set for R 2 . The R 2 values for R2 and R3 were 0.649 and 0.635, respectively. Therefore, no strong dependency on OLR existed. The fits should be described by more complex models that take into account mass transfer relationships. The MY was satisfactorily described by its inverse relationship with OLR for R2. The R 2 was 0.881 with a D-W of 3.6 (p-value = 0.990). While the fit for R3 was smaller, the R 2 of 0.782 with D-W of 2.22 (p-value = 0.4709) was still significant. The inverse relationship between MY and OLR has been described by several authors [44][45][46]. Since both the p-values above were greater than 0.05, there was no indication of serial autocorrelation in the residuals at the 95.0% confidence level. The fitting for both reactors is shown in Figure 7. All the points were contained or at least very close to the confidence limits of the prediction lines (green lines). Therefore, with all the above information combined, the models were considered acceptable.

Optimization of a Reactor
The two equation-systems developed for R2 and R3 combined the Stover-Kincannon model and reciprocal model for MY. Hence, the optimal OLR to simultaneously optimize MPR and MY for R2 and R3 are described by Equations [(11) and (12)] and [ (13) and (14) The graphical optimizations are shown in Figure 8. The ordinates represent the % of the MPRmax or the MYmax measured. The call-out represents the point where both functions met each other. R2 can handle a higher OLR; however, the yields were less from both The graphical optimizations are shown in Figure 8. The ordinates represent the % of the MPR max or the MY max measured. The call-out represents the point where both functions met each other. R2 can handle a higher OLR; however, the yields were less from both functions than R3. Nevertheless, both reactors have similar working intervals, which provided reasonable yields from 3 to 5 gCOD/L for both response variables.
Using the averaged value of the COD in the feed, the working intervals were between 4 and 7 d and 6.5 and 11 d for mixtures 2 and 3, respectively. The upper value for mixture 3 was slightly above 10 d (HRT), the interval for agro-industrial wastewaters suggested by Cruz-Salomón et al. [25]. Meanwhile, mixture 2 had a working interval that was slightly lower than the selected interval. This was attributed to the lower COD content in the mixture. Appl. Sci. 2022, 12, x FOR PEER REVIEW 18 of 23 functions than R3. Nevertheless, both reactors have similar working intervals, which provided reasonable yields from 3 to 5 gCOD/L for both response variables.
Using the averaged value of the COD in the feed, the working intervals were between 4 and 7 d and 6.5 and 11 d for mixtures 2 and 3, respectively. The upper value for mixture 3 was slightly above 10 d (HRT), the interval for agro-industrial wastewaters suggested by Cruz-Salomón et al. [25]. Meanwhile, mixture 2 had a working interval that was slightly lower than the selected interval. This was attributed to the lower COD content in the mixture.

Discussion
The AD of the PM+CWM, PM+CWM+SWW, and PM+CWM+SBT substrate mixtures in a continuous operation mode using the three different EGSBs reactors yielded three main conclusions regarding the mixtures (see below).
1. The synergistic effects described by the batch model in [24] were also found in the continuous operation. 2. The maximum methane yields in the continuous operation of any mixture of these four substrates were predicted using the batch model and multiplying the BMP∞ by a coefficient between 0.7 and 0.8. 3. The employment of the Stover-Kincannon model showed that all three mixtures had a different kinetical behavior, which could even be noticed among the two triple mixtures.
The synergistic effects due to the addition of a third substrate were most likely related to the C/N values. The high C/N values in the three-substrate mixtures explained the good performances observed; however, a higher ratio was no indication of the better performance of mixture 3 versus 2. Hence, the performance of the co-digestion of these mixtures should not be oversimplified by the C/N values without having taken into account other influential factors. Nevertheless, the most recommended C/N values in the literature were from 0 to 30 [4,13,63]. We note that all our mixtures had a C/N value < 20 (Table 1), which strongly suggested the increased proportion of the carbon-rich-substrates within the mixtures.
The concept of integrating an EGSB reactor in a typical agricultural biogas plant is also of relevance. Compared to a typical agricultural biogas plant, where the representative HRT values are between 50 and 150 d [12,64], high-rate reactors provide an alternative system for the treatment of liquid substrates or their liquid fractions at much smaller HRTs. Substrates mixtures that have influenced HRTs should be applied as suggested by Paulose and Kaparaju [63]. They stated that a degradation rate follows an inverse function with the HRT depending on substrate complexity. Consequently, higher HRTs need to be applied and lower degradation rates were expected for lignin-rich substrates than for protein-or sugar-rich substrates. Agricultural biogas plants in Germany typically co-digest animal manure with either maize or grass silage [11]; therefore, higher HRTs are expected for the three-substrate mixtures digested in this paper due to their complexity. However,

Discussion
The AD of the PM+CWM, PM+CWM+SWW, and PM+CWM+SBT substrate mixtures in a continuous operation mode using the three different EGSBs reactors yielded three main conclusions regarding the mixtures (see below).

1.
The synergistic effects described by the batch model in [24] were also found in the continuous operation.

2.
The maximum methane yields in the continuous operation of any mixture of these four substrates were predicted using the batch model and multiplying the BMP ∞ by a coefficient between 0.7 and 0.8.

3.
The employment of the Stover-Kincannon model showed that all three mixtures had a different kinetical behavior, which could even be noticed among the two triple mixtures.
The synergistic effects due to the addition of a third substrate were most likely related to the C/N values. The high C/N values in the three-substrate mixtures explained the good performances observed; however, a higher ratio was no indication of the better performance of mixture 3 versus 2. Hence, the performance of the co-digestion of these mixtures should not be oversimplified by the C/N values without having taken into account other influential factors. Nevertheless, the most recommended C/N values in the literature were from 0 to 30 [4,13,63]. We note that all our mixtures had a C/N value < 20 (Table 1), which strongly suggested the increased proportion of the carbon-rich-substrates within the mixtures.
The concept of integrating an EGSB reactor in a typical agricultural biogas plant is also of relevance. Compared to a typical agricultural biogas plant, where the representative HRT values are between 50 and 150 d [12,64], high-rate reactors provide an alternative system for the treatment of liquid substrates or their liquid fractions at much smaller HRTs. Substrates mixtures that have influenced HRTs should be applied as suggested by Paulose and Kaparaju [63]. They stated that a degradation rate follows an inverse function with the HRT depending on substrate complexity. Consequently, higher HRTs need to be applied and lower degradation rates were expected for lignin-rich substrates than for protein-or sugar-rich substrates. Agricultural biogas plants in Germany typically co-digest animal manure with either maize or grass silage [11]; therefore, higher HRTs are expected for the three-substrate mixtures digested in this paper due to their complexity. However, the differences in the HRTs were always noticeably large. Ruile et al. [12] studied 21 full-scale plants in the region of Baden-Württemberg (southern Germany), which performed either single digestion or co-digestion of cattle manure, maize silage, and grass silage at different solid contents. They found that high values of degradability were reached at a HRT of ≥100 d. Thus, a more sophisticated scheme of treatment that involves multistage processes has been suggested for more efficient energy production [4,63]. Thus, the integration of a high-rate reactor in a typical treatment plant could lead to an increment in energy production, as found by Shen et al. [65] in their co-digestion of fruit/vegetable and food wastes in two stages (UASB+CSTR). This approach allowed them to work at higher OLRs and increase MPR values up to 15% over a single-stage digestion (UASB).
We were able to operate all three reactors for up to 1 d, where the tanks were refilled daily. Consequently, the daily preparation of the mixtures was a logistical and practical challenge. Also, the MYs obtained at a HRT of 1 d were the lowest among the three reactors. It was probably not ideal to run the reactors at such a low HRT; however, this was possible and can be especially useful when the demand for biogas is peaking or excess amounts of substrates need to be processed.
The results of the continuous operation were significantly influenced by the lack of proper mixing in the storage tanks and the seasonal behavior of the substrates. The latter was more obvious in the manure substrates and buffered by the addition of a third substrate. Mixture 2 was the least affected since it had the lowest ODM content in the feeds. Consequently, this mixture had the highest substrate homogeneity inside the reactor due to mixing by recirculation and increased biomass-substrate contact, which facilitated the operation [19,20].
The obtained results support the technical feasibility of the AcoD of liquid manurebased mixtures using EGSB reactors. Thus, it raises the possibility of designing new treatment concepts employing EGSB reactors for the AD of liquid agro-industrial mixtures while significantly reducing the required operating time.
The calculation of the hydraulic dimensionless numbers strongly suggested a CSRT behavior. The assumption of a single reactor with a CSTR behavior simplified the modeling of an EGSB reactor. This was strongly considered when applicable. Due to the lack of biomass sampling along the reactor, and without an adequate computer flow dynamics (CFD) model or RTD study, the consideration of one CSTR seemed the better option. However, the measurements of biomass concentration together with CFD modeling or an RTD study were highly recommended to thoroughly model a reactor [3,6].
Likewise, the operation intervals for an operation at a pilot plant scale were laid down for the two three-substrate mixtures, since both mixtures were likely more profitable than the two-substrate mixture considered. This was in the context of EEG. Furthermore, we recommend the development of more complex models which will allow for the simultaneous control of several process variables as well as describe the potential interactions involved within these variables.

Conclusions
The AcoD of the liquid fraction of PM+CWM and a third carbon-rich substrate such as SWW or SBT was successfully carried out in EGSB reactors, which were operated continuously. This work provides an alternative to typical CSTR systems used for manure and liquid manure treatments. The flow pattern of the studied reactors behaved similarly to a complete mixture reactor. Notably, the hydraulic behavior of our reactors was similar to those found in the literature. Moreover, the results from the batch test were successfully transferred to a continuous scale through the development of empirical and statistical modeling and the optimization of operating OLR intervals. We will consider more complex mechanistic models in the future. Further experiments are going to be carried out at the pilot plant scale in the Saerbeck bioenergy park using one automatically controlled 500 L EGSB reactor.

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

Conflicts of Interest:
The authors declare no conflict of interest.

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

Conflicts of Interest:
The authors declare no conflict of interest.