Methane Concentration Prediction in Anaerobic Codigestion Using Multiple Linear Regression with Integrated Microbial and Operational Data

1. Introduction

The Valle del Cauca region is one of Colombia’s most active agro-industrial areas, combining high agricultural productivity with unique ecological richness. The territory is sustained by ecosystems that range from coastal plains to montane forests, which support both biological diversity and productive capacity. It ranks as the third-largest producer and consumer of pork in the country, with a reported output of 88105 tons in 2023, equivalent to 15.6% of national production, and an average pig population exceeding 396000 animals [1,2]. This sector generates large volumes of pig manure (PM) that require appropriate handling to prevent environmental and public health risks. Another productive activity with growing regional relevance is cassava cultivation, which covered approximately 564 hectares in 2020, yielding a total of 9888 tons of fresh roots [1]. During starch extraction, each kilogram of cassava generates about 0.2 kilograms of starch, 0.65 kilograms of fibrous residue (cassava dregs (CD)), and between five and seven liters of wastewater [3]. Based on these ratios, the estimated annual generation of by-products in the region reaches nearly 6427 tons, most of which are not currently valorized [4,5,6].

To address the increasing accumulation of organic residues from pig farming and cassava processing, anaerobic digestion (AD) has been promoted in rural areas of Valle del Cauca as a strategy for energy recovery and waste management. In these settings, one phase tubular biodigesters are commonly employed due to their affordable construction, ease of installation and minimal infrastructure requirements, making them particularly attractive to smallholder producers [7,8]. AD is a biologically mediated process in which organic matter is sequentially transformed into methane through four main stages, each driven by specific microbial groups. In the hydrolysis phase, hydrolytic bacteria degrade complex macromolecules such as carbohydrates, proteins and lipids into soluble monomers. These compounds are then metabolized by acidogenic bacteria during acidogenesis, producing volatile fatty acids (VFA), alcohols and gases like hydrogen and carbon dioxide. In the subsequent acetogenesis stage, acetogenic microorganisms convert these intermediates into acetate, which, along with hydrogen and carbon dioxide, is used by methanogenic archaea in the final phase to generate methane [9].

For the process to remain stable and efficient, environmental conditions such as pH and temperature must be kept within optimal ranges, typically between 6.5 and 7.5 for pH and 30 to 38 °C under mesophilic conditions [10,11,12,13]. In addition, maintaining a C:N ratio between 20:1 and 30:1 is considered ideal for AD, as it ensures sufficient nitrogen for microbial growth without leading to ammonia inhibition or carbon limitation [14,15,16,17]. However, most rural systems lack monitoring tools and operate through empirical practices, without clear understanding of internal conditions or microbial dynamics [7]. This limitation frequently leads to process imbalance, reduced performance and early system failure.

To overcome the performance limitations of conventional digesters, several strategies have been developed to improve substrate biodegradability and enhance biogas production. Among them, mechanical pre-treatments, codigestion, and multiphase configurations have proven to be particularly effective in increasing system efficiency [18,19,20,21,22]. Mechanical pre-treatments have proven effective in enhancing the hydrolysis of lignocellulosic substrates by reducing particle size and fiber crystallinity, thus increasing surface area and enzymatic accessibility [23,24]. Depending on specific conditions, methane production improvements of 16% to 99% have been reported with mechanical treatments [25]. These results highlight the potential of simple mechanical treatments to enhance biodegradability and biogas productivity, especially during the hydrolysis and acidogenesis phases, which are often rate-limited in solid waste digestion.

Codigestion has emerged as a robust strategy to address the nutrient imbalances and low biodegradability often associated with single-substrate digestion. By combining complementary feedstocks, this approach improves the carbon to nitrogen (C:N) ratio, dilutes inhibitors, and stimulates microbial activity, allowing for higher energy yields [26]. For instance, it has been reported that co-digesting PM with cassava pulp at inclusion levels of up to 60% of the incoming volatile solids (VS) can increase the specific methane yield by around 41% compared with PM alone [27]. Similarly, results have also shown that mixtures containing 66% PM, 16% cassava pulp, and 16% bagasse achieve higher methane yields than those with high bagasse content alone, which led to pH imbalances and process failure [28]. Likewise, experimental trials combining sewage sludge with food waste reported an increase in methane yield from 159 to 799 mL CH₄/gVS, along with a reduction in hydraulic retention time (HRT) to less than four days [29]. These improvements are attributed to enhanced microbial synergy and substrate availability, which accelerate volatile solids degradation.

Multiphase AD systems have been developed to address the limitations of single-stage configurations by creating distinct operational environments for each metabolic phase. In two-phase systems, the acidogenic and methanogenic stages are physically separated, which enables more efficient substrate conversion, greater resilience to organic shocks, and better pH control [30,31]. This structural decoupling has led to increases in methane yields, improved volatile solids removal, and significant reductions in HRT without compromising performance [31]. Compared to single-stage systems, which often suffer from suboptimal compromises between the needs of different microbial groups, two-phase configurations facilitate the coexistence of specialized communities under more stable conditions [31,32]. Although three-phase systems further refine process compartmentalization by isolating hydrolysis, acidogenesis, and methanogenesis, they often entail higher operational complexity, energy consumption, and maintenance requirements [33]. These drawbacks have limited their scalability, particularly in low-resource contexts. Consequently, two-phase systems strike a practical balance between performance enhancement and technical feasibility, making them a more accessible alternative for decentralized applications.

Among the strategies developed to improve AD performance, the integration of real-time monitoring systems has become increasingly relevant for enhancing process oversight and operational efficiency [34,35]. Basic and key variables such as pH, temperature, and methane concentration can be considered to infer the internal state of the reactor and anticipate potential imbalances. The use of cost-effective IoT platforms such as ESP32 microcontrollers coupled with sensor has proven suitable for real-time tracking, achieving deviations below 2% for CH₄ and 1.7% for pH when compared to laboratory-grade methods [34,35,36]. Systems incorporating the MQ-4 sensor (200–10,000 ppm CH₄) and platforms like ThingSpeak facilitate continuous data acquisition, cloud visualization, and automatic alerts, offering a practical solution to reduce manual intervention and increase system reliability [35,37,38].

In parallel, greater attention should be given to the microbial community (MC) involved in AD, as they are rarely considered in routine operation despite being responsible for driving the entire process [9]. Recent studies have highlighted that variations in microbial structure are strongly influenced by substrate type, operational parameters such as temperature and organic loading rates (OLR), and reactor configuration [39,40]. However, most operational strategies still rely exclusively on physicochemical parameters, overlooking microbial signals that often precede system imbalances [41]. Sequencing platforms like Illumina MiSeq and NextSeq have revealed both dominant and low-abundance taxa with key metabolic roles, including members of Euryarchaeota involved in methanogenesis and syntrophic bacteria mediating VFA conversion [40]. Understanding microbial shifts under stress conditions has provided valuable insights into process behavior and system dynamics, reinforcing the need to integrate microbial data into process understanding, particularly to elucidate how shifts in community structure and function impact methane levels [42,43].

Despite their central role in AD, MLR models have traditionally been developed using operational variables that capture external system conditions, parameters that are directly measurable or predefined during setup, while MC have often been treated as secondary inputs or excluded altogether. For example, recent studies have used MLR to predict specific methane production from dry AD of the organic fraction of municipal solid waste in pilot-scale plug-flow reactors. Six significant, mostly operational predictors were prioritized (VS, OLR, HRT, C/N ratio, lignin content, and VFA) via Pearson correlation and PCA. Simple regression showed low performance (R²= 0.3), while the full MLR reached R²= 0.91. A reduced model with four uncorrelated variables (VS, OLR, C/N ratio, lignin content) maintained strong accuracy (R²= 0.87) with fewer inputs. [44]. Similarly, MLR has been applied to predict VFA concentrations in AD of primary and secondary sludge using operational and physicochemical inputs. The model achieved R² values above 0.85 in several scenarios, offering high interpretability and low computational demand. Although less accurate than leading ensemble methods, MLR remains suitable for applications that require clear interpretation of variable influence [45].

Unlike models based solely on operational parameters, an MLR model using the relative abundances of archaeal and bacterial OTUs predicted methane production rates in 149 anaerobic digesters, explaining 66% of the variance with a standard error of 0.12 LCH₄/L·day. Through the integration of PCA, ANOVA, and cubic smoothing, it captured operational effects such as ammonia toxicity and achieved R² values of 0.93 for biogas generation and 0.95 for COD removal across laboratory, pilot, and industrial scales, underscoring its worth in translating complex MC data into scalable process predictions. [46].

Building on emerging evidence supporting the integration of microbial data into statistical modeling, this study evaluates the potential of MLR to predict methane concentrations in a low-cost, two-phase anaerobic digester treating PM and CD at laboratory scale. The work aligns with Sustainable Development Goal 7 by promoting accessible tools for energy recovery from organic waste.

This article is structured into four main sections. The Introduction outlines the context of AD in the Valle del Cauca region, highlighting environmental and operational challenges from agro-industrial organic waste, reviewing strategies to improve biogas systems, and emphasizing the need to integrate microbial data into predictive models. The Materials and methods detail the system setup, monitoring, sequencing, and the MLR approach used for variable selection and model construction. The Results and discussion section presents the modeling outcomes, identifies relevant predictors, and interprets their contribution to system behavior. The Conclusion summarizes the key findings and future perspectives for incorporating microbiota into data-driven frameworks for sustainable energy transitions.

2. Materials and Methods

This study aims to develop a predictive model for methane concentration based on a set of measurable variables, including VFAs, microbial populations, and operational parameters. This section first describes the dataset and the preprocessing steps undertaken. Subsequently, it details the initial linear modeling approach, followed by a feature selection process based on variable weighting to derive a simplified, yet robust, model. Finally, it presents the development of an adaptive predictive model using a moving window technique combined with a regularization method to prevent overfitting.

2.1. Substrate Selection

The substrates used in this study were fresh PM and CD. The inoculum, obtained from the same source as the manure, was included to ensure microbial compatibility with the feedstock. Both were collected at a small-scale pig farm located in the municipality of Florida, Valle del Cauca, where approximately 20 pigs are kept under semi-intensive conditions. Animal pens are washed twice daily, and the resulting wastewater, rich in organic matter, drains into a static open-air tank that served as the inoculum source. Fresh manure was manually collected after excretion using sanitized tools. CD were obtained from a medium-sized cassava starch-processing facility located in the rural area of Mandiba, Santander de Quilichao, Cauca. Processing nearly eight tons of cassava per day, the plant generates over two tons of lignocellulosic residue each week. This material was delivered in dry, milled form.

All samples were stored at 4 °C until physicochemical characterization, which included proximate analysis by gravimetric methods and determination of the carbon-to-nitrogen (C:N) ratio via high-temperature combustion. These procedures followed the Standard Methods for the Examination of Water and Wastewater (APHA, AWWA, WEF), ensuring analytical consistency as summarized in Table 1 [47].

2.2. Experimental Setup

The experimental setup consisted of a two-phase laboratory-scale anaerobic digester designed to operate without integrated control systems Figure 1. The system was constructed using 110 mm sanitary-grade PVC tubing due to its low cost, durability, and ease of assembly. Phase 1 (D1F1) (3 L) was expected to perform hydrolysis and acidogenesis, while phase 2 (D1F2) (4 L) supposedly supported acetogenesis and methanogenesis. Each chamber was operated at 80% of its total volume, 2.4 L in phase 1 and 3.2 L in phase 2, leaving the remaining headspace for biogas accumulation. To enable real-time monitoring, a low-cost IoT module was incorporated into the digester, integrating an Arduino UNO microcontroller with sensors for pH, temperature, and methane concentration. Data were transmitted through a mobile network to the ThingSpeak platform for remote visualization [38]. This setup allowed continuous monitoring without the need for sophisticated instrumentation.

2.3. Operational Parameter

To establish an active MC, both phases were fed inoculum for five days, until its working volume. The inoculum had a C:N ratio of 10.3 and 2.2% TS. During start-up, the OLR, estimated with a five-day HRT, was 8.37 gVS/L·day. Thereafter, feeding used a 73:27 blend of PM and CD. The daily feed was 35 g fresh PM and 13 g CD, plus 166 g water to achieve 10% TS (214 g/day total). The theoretical C:N ratio was 21.55. With the defined working volumes, HRTs were 12 days for D1F1 and 15 days for D1F2. Corresponding OLRs were 7.7 and 5.7 gVS/L·day. VS inputs were 18.46 g/day (D1F1) and 18.45 g/day (D1F2). Daily manual feeding with graduated containers and isolation valves ensured accurate dosing and anaerobiosis.

The IoT-instrumented digester (D1) enabled incremental, data-driven feed adjustments in both phases (D1F1, D1F2) using real-time pH, temperature, and methane concentration. These signals guided when to lower the OLR and TS and when to apply temporary pH control, moving the reactors toward consistent operating conditions. Five feed formulations were implemented (Table 2). In D1F1, pH was briefly corrected with lime and then NaOH to keep it within 6.5–7.5; by mixture 5, recirculated digestate from D1F2 maintained pH without further chemicals. Mixture 4 used inoculum from an anaerobic digester at a university in Colombia treating food waste. Across mixtures, TS was reduced from 10% to 8–9%, OLR decreased from 12.4 gVS/L·day (inoculum step) to 5–6 gVS/L·day, and the C:N ratio increased in the final mixture due to recirculation while the contributions of PM and CD were reduced.

2.4. Steady-State

Identifying steady-state periods was essential to build a reliable dataset, define representative operating conditions, and guide downstream variable prioritization and modeling. pH, temperature, and methane concentration were monitored continuously for 161 days (24/7). The IoT system logged three readings per minute for each variable and was routinely cross-checked against bench measurements to validate operational reliability.

Data volume was substantial, D1F1 recorded 694110 samples per variable and D1F2 573215. Processing followed six steps: 1) splitting timestamp into date and time; 2) validity filtering (e.g., pH 3–12; 10–45 °C; CH₄ within instrument bounds) with out-of-range values set to blank; 3) multivariate imputation by chained equations (MICE) to preserve temporal continuity [48]; 4) resampling to hourly means (2893 rows in D1F1; 2389 in D1F2) and 5) to daily means (152 and 147, respectively), retaining trends while reducing computational load as shown in Table 3.

Stable windows were then identified via rolling windows using relative standard deviation thresholds (<15%) around moving means for pH, temperature, and methane concentration, with a minimum continuous duration and compliance with predefined operating limits [49]. D1 showed extended steady windows, typically with pH 6.5–7.5, facilitated by high-frequency data and the ability to adjust operating conditions in real time.

2.5. VFA Quantification

Samples were collected every three days in 5 ml Eppendorf tubes and stored at -20°C until analysis. The final selection of samples for analysis was made considering the periods of system stabilization under IoT monitoring and budgetary constraints, prioritizing those most representative of the overall process behavior. Sampling was carried out during the active operation of the digester.

The quantification of VFAs was performed by gas chromatography, following the procedure described in section 5560D of the Standard Methods for the Examination of Water and Wastewater (APHA) [50], in the laboratory of the Department of Chemical Engineering and Analytical Chemistry at the University of Barcelona. Prior to chromatographic analysis, the samples were centrifuged and filtered through 0.45 µm nylon membranes to remove suspended solids. Each analysis vial contained 1 ml of sample, diluted or not depending on the estimated concentration level, along with 0.1 ml of 15% orthophosphoric acid containing a known concentration of 2-ethylbutyric acid (~500 mg/L) as an internal standard. This compound allowed verification of injection consistency and facilitated calibration of the equipment through the ratio of analyte to standard peak areas.

Analyses were carried out on a Shimadzu GC-2010 Plus gas chromatograph with a flame ionization detector, using a DB-FFAP capillary column (Agilent Technologies, 30 m × 0.25 mm × 0.25 µm). The oven temperature program started at 60°C with a two-minute hold, followed by an increase of 20°C/min up to 240°C, maintained for an additional two minutes. The total analysis time was 13 minutes.

The injector (SPL-1) operated at 220°C in split mode, with a split ratio of 50:1. Helium was used as the carrier gas at a pressure of 42.6 kPa, with a total flow of 233.4 ml/min, a column flow of 8.86 ml/min, and a linear velocity of 60 cm/s. The purge flow was set at 3 ml/min, and the makeup gas flow (nitrogen) at the detector was 10 ml/min. The injection volume was 2 ml, using helium, air, hydrogen, and nitrogen as auxiliary gases.

For equipment calibration, a commercial VFA standard (Supelco CRM46975) containing defined concentrations of acetic, propionic, isobutyric, butyric, isovaleric, valeric, isocaproic, caproic, hexanoic, and heptanoic acids was used. Serial dilutions were prepared in 1:1, 1:2, 1:4, 1:8, 1:16, and 1:32 ratios, to which orthophosphoric acid and the internal standard were also added. For alcohol analysis (ethanol, propanol, and butanol), defined-concentration standard solutions were prepared, applying the same dilutions and analytical conditions.

This procedure allowed precise and reproducible determination of VFAs in the samples, essential for evaluating the performance of the AD system and its relationship with operating conditions and microbiota.

2.6. Metagenomic Analysis

Samples for metagenomic analysis were collected directly from operational biodigester using 50 ml Falcon tubes. Sampling was performed every three days throughout the process, following the same prioritization criteria used for the quantification of VFAs, focusing on periods of greatest microbiological representativeness and considering the availability of resources. Once collected, samples were immediately frozen at -20°C and stored until further processing.

To analyse the MC, Falcon tubes were sent to Omega Bioservices (U.S.A) for DNA extraction using the kit E.Z.N.A.® Universal Pathogen Kit, library preparation and for sequencing the V3–V4 hypervariable region of the 16S rRNA gene using the primers 341F (CCTACGGGNGGCWGCAG) and 806R (GACTACHVGGGTATCTAATCC) which was conducted on an Illumina Miseq sequencing platform (Paired-end sequencing 300 bp). Illumina reads were then analysed using BaseSpace app (version 1.1.3) [51]. Thus, raw sequence data were demultiplexed and then quality filtered, denoised, merged, and chimera removed using the DADA2 [52] to generate amplicon sequence variants (ASVs). Taxonomic assignment was conducted using the SILVA database (version 138.2) [53].

To structure the analysis of microbial interactions, a subset of phyla of interest was defined from the general metagenomic dataset, considering the sequencing reads obtained for each taxonomic group. The selection was based on two main criteria. First, the sustained presence of each phylum throughout the monitoring period was evaluated, excluding those with very low or intermittent representation, as their variability would hinder the detection of consistent associations in the relational analysis. Second, functional relevance reported in previous studies on anaerobic digestion was reviewed, prioritizing phyla whose involvement in fermentative, acetogenic, or methanogenic pathways has been extensively documented in similar systems [9,54].

Once the representative periods were defined, the results from VFA quantification and metagenomic analysis were integrated, extending the characterization to the biochemical and microbiological components of the system. In several cases, the observed patterns were consistent with those reported in the specialized literature, which supported the robustness of the approach. The dataset included operational, biochemical, and microbiological variables [55,56,57,58,59].

Since the biochemical and microbiological measurements were less frequent than the operational records, imputation techniques were applied within the selected periods to expand the dataset without distorting the relationships among variables. Methods such as KNN imputation, iterative imputation, and MICE were employed [56,60].

The analysis focused on the period between days 97 and 154, which, although not representing a fully stabilized phase, shows a trend toward stabilization and coincides with the selected VFA and microbiological samples. This ensured consistency between the experimental data and the operational conditions.

2.7. Preprocessing and Unified Database

Once the representative periods were defined, the results from VFA quantification and metagenomic analysis were incorporated to extend the characterization of the system to its biochemical and microbiological dimensions. The patterns obtained aligned with those reported in specialized literature, reinforcing the validity of the approach [55,56,57,58,59]. The unified dataset combined operational, biochemical, and microbiological variables.

Because biochemical and microbiological measurements were less frequent than operational records, imputation method MICE was applied to harmonize the dataset without altering the underlying relationships among variables [56,60].

The analysis focused on the period between days 97 and 154, which, while not fully stabilized, displayed a clear trend toward steady performance and coincided with the VFA and microbiological samples selected. This ensured coherence between experimental observations and operational conditions. The resulting dataset comprised daily averages over 58 days, which were further refined through linear interpolation to increase temporal resolution. This process expanded the series to 1000 points, enabling the application of moving window analyses, as illustrated in Figure 2.

The interpolation was validated for all variables, yielding R² values close to 1 and MRE values around 0.1%, confirming a high-fidelity representation of the original data.

2.8. Linear Modeling

To simplify the proposed equations and procedure, the suffixes associated with each fatty acid (Table 4) microorganism (Table 5), and operating condition (Table 6) are shown below.

The equation (1) that linearly approximates $\mathrm{C}{\mathrm{H}}_4$ concentration as a function of the microorganisms, fatty acids, and operating conditions was proposed in the following linear form, based on the suffixes from Table 4, Table 5 and Table 6.

$${\left(C{H}_4\right)}_{aprox}=C_1\left(a_1\right)+C_2\left(a_2\right)+\dots +C_{25}\left(m_{25}\right)+C_{26}\left(m_{26}\right)+\dots +C_{30}\left(p_{30}\right)+C_{31}\left(p_{31}\right)$$

(1)

In its matrix form (matrix $\underset{\_}{A}$), equation (1) can be expressed as follows (equation (2)):

$$\left[\begin{array}{c}{\left(C{H}_4\right)}_1\\ {\left(C{H}_4\right)}_2\\ \begin{array}{c}.\\ .\\ .\\ {\left(C{H}_4\right)}_n\end{array}\end{array}\right]= \left[\begin{array}{ccc}{\left(a_1\right)}_1& {\left(a_2\right)}_1& \begin{array}{cc}\dots & \begin{array}{ccc}{\left(m_{26}\right)}_1& \dots & {\left(p_{31}\right)}_1\end{array}\end{array}\\ {\left(a_1\right)}_2& {\left(a_2\right)}_2& \begin{array}{cc}\dots & \begin{array}{ccc}{\left(m_{26}\right)}_2& \dots & {\left(p_{31}\right)}_2\end{array}\end{array}\\ \begin{array}{c}.\\ .\\ .\\ {\left(a_1\right)}_n\end{array}& \begin{array}{c}.\\ .\\ .\\ {\left(a\right)}_n\end{array}& \begin{array}{c}.\\ .\\ .\\ \begin{array}{cc}\dots & \begin{array}{ccc}{\left(m_{26}\right)}_n& \dots & {\left(p_{31}\right)}_n\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ \begin{array}{c}.\\ .\\ .\\ C_{31}\end{array}\end{array}\right]$$

(2)

where the constants $C_i$ are the approximation coefficients. This matrix form from equation (2) can be written more compactly as shown in equation (3):

$$\underset{\_}{C{H}_4}= \underset{\_}{A} \underset{\_}{x}$$

(3)

The matrix $\underset{\_}{A}$ contains data collected from fatty acids, microorganisms, and operating conditions, the vector $\underset{\_}{C{H}_4}$ represents the collected methane production data, while the vector $\underset{\_}{x}$ contains the approximation coefficients that must be determined to formulate the model. The vector $\underset{\_}{x}$, can be solved by rearranging equation (3) as follows:

$$\underset{\_}{x}= [ {\underset{\_}{A}}^T \underset{\_}{A} {]}^{-1} {\underset{\_}{A}}^T \underset{\_}{C{H}_4}$$

(4)

where ${\underset{\_}{A}}^T$is the transpose of matrix $\underset{\_}{A}$.

2.8.1. Assessing Variable Importance

To determine the relative importance of each variable in the approximation, and subsequently define a smaller, more practical subset (as working with all 31 variables can be impractical and costly in terms of laboratory testing), a variable weighting method was used. Therefore, equation (5) appears as a modification of the equation (3) considering the minimal error $ϵ$.

$$\underset{\_}{C{H}_4}=\underset{\_}{A} \underset{\_}{x}+ϵ$$

(5)

To quantify how much each variable "contributes" to the $C{H}_4$ production within the approximation, it is necessary to measure the relevance of each variable in the linear model. Since each variable may be measured on a different scale (e.g., microorganism abundance vs. fatty acid concentration in mg/L), directly comparing the raw coefficients $C_i$ in vector $\underset{\_}{x}$ can be misleading. Therefore, it is necessary to standardize the input data. In the same way, to compare the relative importance of each variable, the coefficients $C_i$ that form the vector $\underset{\_}{x}$ were standardized (as z-scores). The standardized coefficient $C_i^{*}$ for each variable $C_i$ was calculated as:

$$C_i^{*}=C_i{\displaystyle \frac{{\sigma }_{C_i}}{{\sigma }_{CH4}}}$$

(6)

where ${\sigma }_{C_i}$ and ${\sigma }_{CH4}$ are the standard deviations of the approximation coefficient $C_i$ and the response variable $C{H}_4$, respectively.

2.8.2. Proposing a Predictive Model

To capture the evolutionary nature of the anaerobic digestion process, a dynamic predictive model was developed based on a moving windows approach. The model operates iteratively. At each time step $t$, a linear regression model is trained using a window containing the last $k$ observations (in this case, $k=10$ was chosen). This model is then used to make a one-step-ahead prediction of $\underset{\_}{C{H}_4}$ (denoted as $\underset{\_}{\widehat{C{H}_4}}$), as a function of the weighted variables previously described.

However, the use of small data windows can lead to overfitting. To address this problem and improve the model's generalization capability, Ridge Regression was used instead of ordinary least squares. This regression introduces a penalty term into the least squares cost function. For each time window, the objective is to find the coefficient vector $\underset{\_}{x}$ that minimizes the following function:

$${min}_{\underset{\_}{x}}\left({‖\underset{\_}{C{H}_4}- \underset{\_}{A} \underset{\_}{x}‖}_2^2+\lambda {‖\underset{\_}{x}‖}_2^2\right)$$

(7)

where ${‖\underset{\_}{C{H}_4}- \underset{\_}{A} \underset{\_}{x}‖}_2^2$ is the sum of squared errors (the data fit term at time $t$), $\lambda {‖\underset{\_}{x}‖}_2^2$ is the regularization term applied at time $t$, $\lambda$ is the regularization hyperparameter that controls the balance between the data fit and the model simplicity, and $\underset{\_}{x}$ contains the coefficients $C_i$. The hyperparameter $\lambda$ was selected to improve the model predictive performance ($\lambda =2$). Thus, equation (4) is rewritten to obtain the predictive parameters ($\underset{\_}{\widehat{x}}$) by solving the following equation:

$${\left(\underset{\_}{\widehat{x}}\right)}_{t+1}={\left([ {\underset{\_}{A}}^T \underset{\_}{A} + \lambda \underset{\_}{I} {]}^{-1} {\underset{\_}{A}}^T \underset{\_}{C{H}_4}\right)}_t$$

(8)

where $\underset{\_}{I}$ is the identity matrix. The goal of equation (8) is to find the value of the coefficients in vector $\underset{\_}{x}$ at time $t+1$ using the data available at time $t$. Using equation (8), it is possible to find the $\underset{\_}{\widehat{C{H}_4}}$ values for a subsequent window, given a defined window size of $k=10$. In this way, equation (3) becomes a prediction equation as follows:

$${\left(\widehat{\underset{\_}{C{H}_4}}\right)}_{t+1}= {\underset{\_}{\left(A\right)}}_{t+1} ( {\widehat{\underset{\_}{x}} )}_{t+1}$$

(9)

2.8.3. Model Performance Evaluation

The precision of the predictive model was quantified using three standard statistical metrics. These metrics evaluate the divergence between the real observed values of $\underset{\_}{C{H}_4}$ and the values predicted by the model, $\widehat{\underset{\_}{C{H}_4}}$. On one hand, the Coefficient of Determination ($R^2$) indicates the proportion of the variance in methane production that is predictable from the independent variables. A value close to 1 indicates an almost perfect fit. Equation (10) shows how it was calculated for this case.

$$R^2=1-{\displaystyle \frac{\sum _{i=1}^n{\left[{\left(C{H}_4\right)}_i-{\left(\widehat{C{H}_4}\right)}_i\right]}^2}{\sum _{i=1}^n{\left[\left({\left(C{H}_4\right)}_i-\overline{C{H}_4}\right)\right]}^2}}$$

(10)

Next, the Root Mean Square Error (RMSE) represents the standard deviation of the prediction residuals. It is a measure of the average error of the model in the same units as the response variable (ppm of $C{H}_4$), which facilitates its interpretation and is expressed in equation (11).

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}} \sum _{i=1}^n{\left[{\left(C{H}_4\right)}_i-{\left(\widehat{C{H}_4}\right)}_i\right]}^2}$$

(11)

Finally, the Mean Relative Error (MRE) measures the average error in relative or percentage terms with respect to the real value is defined in equation (12). The absolute value was used to prevent positive and negative errors from canceling each other out.

$$MRE={\displaystyle \frac{1}{n}} \sum _{i=1}^n\left|{\displaystyle \frac{{\left(C{H}_4\right)}_i-{\left(\widehat{C{H}_4}\right)}_i}{{\left(C{H}_4\right)}_i}}\right|$$

(12)

In equations (10), (11), and (12), ${\left(C{H}_4\right)}_i$ is the real value of the i-th observation, ${\left(\widehat{C{H}_4}\right)}_i$ is the value predicted by the model for the i-th observation, $\overline{C{H}_4}$ is the mean value of all real values, and n is the total number of observations used for the evaluation.

3. Results and Discussion

3.1. Digester Performance

During the implementation of the two-phase biodigester, one of the main technical challenges was controlling gas leaks and internal pressure, which required multiple structural adjustments and caused delays in the early stages of operation. The installed manometers failed, likely due to H₂S-induced corrosion, while the valves progressively stiffened with use, in some cases requiring replacement. In addition, certain PVC weld joints developed fractures, compromising system integrity. The manual stirring mechanism did not improve process performance but was associated with gas leaks, leading to its deactivation and reinforcement of seals. After several corrective interventions, continuous and functional operation was achieved for the duration of the experiment.

3.1.1. IoT Monitoring Advantages

Although low-cost systems might be unsuited for long-term use or high-precision data collection, these biodigesters provide a practical alternative for experimental applications at laboratory scale when resources are limited. The estimated cost of assembling a two-phase biodigester without IoT monitoring was 300 USD, while the addition of a digital monitoring system increased the total to 420 USD per unit. Although some sensors required replacement during operation, the IoT system performed reliably, providing consistent readings comparable to manual instruments. Its implementation enabled real-time, continuous data acquisition, which was essential for detecting operational variations, making timely adjustments, and improving process understanding.

Figure 3 (D1F1) and Figure 4 (D1F2) show the 24-hour profiles of pH, temperature, and CH₄ concentration during different operational stages. The days were randomly selected from both phases to provide representative snapshots of system behaviour under varying conditions. In all cases a consistent inverse relationship was observed, where pH increased during early morning hours as ambient temperature declined and then decreased progressively as temperature rose throughout the day. In D1F1 this trend was stronger and more reproducible, with correlation coefficients between -0.84 and -0.94, while in D1F2 the association was weaker (between -0.49 and -0.66) and accompanied by larger fluctuations in methane concentrations. These results highlight the direct effect of ambient thermal oscillations on microbial activity, especially on pH dynamics.

This behaviour may be linked to phases of microbial adaptation or to the accumulation of internal self-regulation mechanisms. As temperature dropped during the night and metabolic activity slowed, nitrogenous compounds likely continued decomposing and releasing ammonia (NH₃). This ammonia could react with dissolved CO₂, which is more soluble at low temperatures, to form ammonium bicarbonate (NH₄HCO₃). The resulting increase in alkalinity buffered pH variations, preventing excessive acidification and contributing to system resilience [17,61].

Such dynamics are rarely captured in conventional laboratory-scale digesters where data are typically restricted to discrete measurements. In this case the use of IoT-based continuous monitoring provided hourly resolution and made it possible to identify fine-scale responses such as the rise in pH at lower nighttime temperatures that would otherwise remain unnoticed. This approach delivers a more realistic picture of system performance under environmental conditions and emphasizes the value of continuous monitoring strategies for interpreting anaerobic digestion behaviour beyond the limits of punctual sampling.

3.1.2. Stabilization of Anaerobic Codigestion

Achieving a steady-state is a critical milestone in AD, as it reflects the convergence of operational and microbial conditions that support sustained methanogenic activity [62,63]. In system D1, specific time segments were identified where pH, temperature, and methane concentration aligned within the functional ranges expected for AD. Figure 5 integrates these three variables across the full operational period, providing a comprehensive view of the transitions from unstable to stabilized phases. This visualization not only highlights the progression of the system under different operating conditions but also illustrates how corrective measures and phase-specific dynamics gradually steered the reactor toward a functional equilibrium.

In Figure 5a, a shift becomes evident after day 120, when pH consistently remained above 6.5 while methane rose steadily, surpassing 8000 ppm by day 132. These conditions coincided with stable temperatures between 30 and 32 °C, an optimal mesophilic range that favours methanogenic activity [49,55]. The segmentation into D1F1 and D1F2 reveals the influence of each phase under a combined HRT of 27 days. During early mixtures (Mix1-Mix2), high OLR and the absence of adapted inoculum produced irregular methane signals dominated by acidogenesis. From Mix3 onwards, corrective measures such as alkalinization promoted higher CH₄ concentrations, although fluctuations beyond ±15% prevented these periods from being classified as steady. Toward Mix4 and Mix5, adjustments including higher inoculum input and recirculation from D1F2 likely increased microbial density and functional diversity, progressively creating conditions more favourable to methanogenesis [64,65].

Figure 5b highlights two segments where the system approached steady-state behaviour. The first, between days 126 and 136, was characterized by pH values between 6.5 and 7.5, methane above 8000 ppm, and stable temperatures (30-32 °C), all within functional ranges and with fluctuations below ±15%. A temporary decline in methane around day 137 disrupted this stage, but from day 139 onwards the system recovered, initiating a second steady segment that persisted until the end of the experiment.

Overall, the convergence of pH, CH₄, and temperature, demonstrates that the steady-state achieved in D1 was not the result of a single correction but the cumulative effect of progressive adjustments. This sequence of changes allowed the system to transition from acidogenic predominance to a consolidated methanogenic phase, representing a functional stabilization consistent with the goals of two-phase anaerobic digestion [65].

3.2. Volatile Fatty Acids (VFAs) and Metagenomic Analysis

Once stabilization was established from IoT-monitored variables, VFAs and microbiota were analyzed during the transition toward optimal operation (days 97-154, Mix 4 and Mix 5). Thirteen samples were taken, and missing data were inputted through MICE to ensure continuity.

Figure 6 shows that cumulative VFA concentrations dropped sharply after day 100, from above 14000 mg/L to 5800 mg/L, before oscillating between 5000 and 7000 mg/L. This decline reflects the mitigation of acidogenic pressure and the progressive adjustment of the microbial community, setting conditions increasingly suitable for methanogenesis and linking metabolite dynamics with microbial responses in the path toward functional balance (Figure 8) [66,67].

The individual analysis of VFAs confirmed that the steep decline after day 100 was largely driven by the reduction of acetic and butyric acids, both tied to early fermentative pathways [66]. Between days 103 and 118, however, propionic acid and medium-chain carboxylates (C5-C8), including caproic, heptanoic, and valeric, increased notably, reaching averages of 140 mg/L, 218 mg/L, and 820 mg/L, respectively [68]. These less common metabolites are typically linked to secondary fermentation processes or to transitional phases of temporary accumulation [69,70,71]. Their persistence, together with measurable levels of propanol (134 mg/L) and the absence of ethanol, suggests a fermentative stage dominated by chain-elongation routes, potentially hindered by propanol’s inhibitory effect on methanogenic consortia [69,72,73]. After day 119, these acids gradually declined (e.g., caproic down to 109 mg/L, valeric to 758 mg/L), while ethanol reappeared (20 mg/L) and propanol rose to 191 mg/L. This pattern may indicate that, despite higher alcohol concentrations, microbes capable of degrading medium-chain acids regained activity, backing a functional shift toward methanogenesis [74,75].

Regarding microbial analysis, a total of 1815465 high-quality reads, with an average of 201718 ± 103945 reads per sample. Rarefaction analysis based on Shannon index showed that sequencing depth was adequate to capture most of the bacterial diversity across samples as shown in Figure 7.

The microbial dynamics based on the phylum level aggregate counts (Supplementary Table S1) reflected in Figure 8 provide insights into the links between dominant phyla and the VFA and methane concentration profiles presented in Figure 5 and Figure 6. During days 97 and 154, Firmicutes remained the prevailing group, averaging 49235 reads (64.2% of the total), underscoring its central role in the early stages of the process, particularly in hydrolysis and acidogenesis [9]. This activity likely promoted the production of fermentative precursors, consistent with the elevated provide insights into the links between dominant acids recorded at the beginning of this interval [66].

Along with Firmicutes, Bacteroidetes (15%) and Actinobacteria (7.5%) contributed to medium-chain fatty acids such as valeric and caproic during the accumulation phase (days 103-118) [69,73]. This functional diversity points to a bacterial consortium engaged in degrading complex polymers and extending fermentative pathways, buffering intermediates before methanogenic activity resumed [76,77]. Toward the end, Firmicutes declined while Euryarchaeota increased to 1.5%, coinciding with reduced VFAs and steadier methane, suggesting activation of acetoclastic and hydrogenotrophic routes [9,78].

Minor groups might have played complementary roles, with Planctomycetes (3.9%) likely coupling sulfide oxidation to methanogenesis, Proteobacteria (1.9%) contributing to propionate and acetate turnover, and Synergistetes (1.9%) participating in syntrophic H₂ transfer [9,79].

Unexpectedly, Verrucomicrobia appeared in the community profile, a phylum typically restricted to volcanic habitats dominated by acidophilic methanotrophs [80]. These bacteria can oxidize methane as their main substrate and, to a lesser extent, hydrogen, carbon dioxide, ammonium and hydrogen sulfide, functioning as natural biofilters in extreme ecosystems. Their presence in a mesophilic anaerobic digester is unusual and may reflect residual inoculum or localized microredox niches rather than an active role in methanogenesis [80,81]. Alongside other low-abundance phyla such as Lentisphaerae, Candidatus saccharibacteria, and Parcubacteria, their detection expands the taxonomic spectrum and raises questions about potential ecological roles still unexplored in anaerobic bioenergy systems [82]. Many of these groups remain unresolved at the species level, even after advanced genomic assembly, forming part of the so called microbial dark matter. This hidden fraction highlights one of the major challenges in deciphering the functional complexity of anaerobic microbiomes [83,84].

3.3. Multiple Linear Regression (MLR)

Modelling phase was based on the Supplementary Table S2. he values in Table 7 show the coefficients ${\mathrm{C}}_{\mathrm{i}}$obtained by finding the vector $\underset{\_}{\mathrm{x}}$ after applying equation (4).

This method of finding the approximation coefficients is known by some authors as inverse modeling and can be considered a linear regression. When the vector $\underset{\_}{x}$ was obtained, equation (1) was applied to generate the approximation curve. Figure 9 shows the resulting approximation.

3.3.1. Data Prioritization

While Figure 9 shows the overall approximation result, it is important to determine the contribution of each fatty acid, microorganism, or operating condition to the $C{H}_4$ production. As previously mentioned, a direct comparison of the $C_i$coefficients can be misleading, so it is important to perform a variable weighting process using equation (6). The values obtained from this process are shown in Table 8 and plotted in Figure 10.

To define how many variables are needed to achieve a good fit without a significant loss of precision, two limits were established, as can be seen in Figure 11 : an ${\mathrm{R}}^2$ Coefficient greater than 0.9 and an MRE lower than 15%. Based on Table 8 and Figure 10, the 12 variables with the highest values, or greatest impact on $C{H}_4$ production, were selected. There are $C_{10}^{*}$, $C_{12}^{*}$, $C_{17}^{*}$, $C_{13}^{*}$, $C_{15}^{*}$, $C_{11}^{*}$, $C_{20}^{*}$, $C_{14}^{*}$, $C_9^{*}$, $C_{16}^{*}$, $C_{19}^{*}$ and $C_{18}^{*}$. According to Table 5, these correspond to Bacteroidetes, Proteobacteria, Verrucomicrobia, Planctomycetes, Spirochaetes, Actinobacteria, Armatimonadetes, Synergistetes, Firmicutes, Euryarchaeota, Tenericutes, and Cloacimonetes, respectively.

Thus, equations (1) and (2) can be rewritten in terms of the 12 variables that were found to be most important according to the weighting performed. This leads to equation (13):

$$\begin{gathered} {\left(C{H}_4\right)}_{WeightedAprox}=C_{10}\left(m_{10}\right)+C_{12}\left(m_{12}\right)+C_{17}\left(m_{17}\right)+C_{13}\left(m_{13}\right)+ \\ {C}_{15}\left(m_{15}\right)+ C_{11}\left(m_{11}\right)+ C_{20}\left(m_{20}\right)+ C_{14}\left(m_{14}\right)+C_9\left(m_9\right)+ C_{16}\left(m_{16}\right)+ \\ {C}_{19}\left(m_{19}\right)+C_{18}\left(m_{18}\right) \end{gathered}$$

(13)

In its matrix form, equation (13) can be expressed as follows:

$$\left[\begin{array}{c}{\left(C{H}_4\right)}_1\\ {\left(C{H}_4\right)}_2\\ \begin{array}{c}.\\ .\\ .\\ {\left(C{H}_4\right)}_n\end{array}\end{array}\right]= \left[\begin{array}{ccc}{\left(m_{10}\right)}_1& {\left(m_{12}\right)}_1& \begin{array}{cc}\dots & {\left(m_{18}\right)}_1\end{array}\\ {\left(m_{10}\right)}_2& {\left(m_{12}\right)}_2& \begin{array}{cc}\dots & {\left(m_{18}\right)}_2\end{array}\\ \begin{array}{c}.\\ .\\ .\\ {\left(m_{10}\right)}_n\end{array}& \begin{array}{c}.\\ .\\ .\\ {\left(m_{12}\right)}_n\end{array}& \begin{array}{c}.\\ .\\ .\\ \begin{array}{cc}\dots & {\left(m_{18}\right)}_n\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{C}_{10}\\ C_{12}\\ \begin{array}{c}.\\ .\\ .\\ C_{18}\end{array}\end{array}\right]$$

(14)

Table 9 shows the new values found for the approximation coefficients, calculated using equation (13) in conjunction with equation (4). These results in the weighted approximation shown in Figure 12.

Table 10 shows a comparison of the metrics used to compare the fits discussed previously. While the general behavior is replicated by both curves, the $R^2$ coefficient for the approximation curve using equation (1) is 0.989, whereas with the weighted approximation from equation (13), the $R^2$ value is 0.979.

Meanwhile, the mean relative error (MRE) for the approximation with equation (1) is around 12.59%, while for the weighted approximation with equation (13), it is 14.94%. Finally, the errors evaluated by the RMSE are below 450 ppm, which, considering the scale of Figure 5, are within an acceptable range. This approach successfully developed a simplified, dynamic model for predicting methane ($C{H}_4$) production in an anaerobic digestion process. The key achievement was the ability to reduce a complex system of 31 variables to a robust predictive model based on only the 12 most influential factors, without a significant loss of precision.

The primary finding of this work is the overwhelming importance of microbial populations as indicators of $C{H}_4$ concentration compared to VFAs and operational parameters. The variable weighting analysis revealed that the 12 most significant variables were exclusively microorganisms, with groups like Bacteroidetes and Proteobacteria showing the highest importance scores. This suggests that, within the context of this study, the state of the microbial community is a more direct and powerful predictor of methanogenic activity than the concentration of intermediate substrates (VFAs) or the operational conditions measured. While VFAs are essential for methanogenesis, their concentrations can be transient. In contrast, the abundance of specific microbial groups likely represents the metabolic potential of the system, making them more robust indicators for modeling purposes.

3.3.2. Predictive Model

The simplification of the model from 31 variables to 12 demonstrates the practical value of the feature selection process. The weighted approximation model, using only the selected microorganisms, achieved an $R^2$ of 0.979, a negligible decrease from the 0.989 $R^2$ of the full model. By focusing only on the most critical microbial indicators, laboratory testing and data analysis efforts can be substantially reduced while still maintaining a high degree of accuracy. The slight increase in MRE and RMSE is an acceptable trade-off for the considerable reduction in model complexity.

In this context, using equation (13), it was possible to develop a predictive model to predict the components of the vector $\underset{\_}{\widehat{x}}$ and, subsequently, the behavior of the $\underset{\_}{\widehat{C{H}_4}}$ production by applying equations (8) and (9), respectively. This was done by using 70% of the dataset to train the model and the remaining 30% to test its performance. Thus, Figure 13 shows how, by applying equation (8), it is possible to obtain a prediction for the behavior of each of the coefficients for the weighted variables $\underset{\_}{\widehat{x}}$ (originally listed in Table 9).

By applying equation (9), the prediction for $\widehat{\mathrm{C}{\mathrm{H}}_4}$ was obtained, as shown in Figure 14. This was derived from the coefficients $\underset{\_}{\widehat{x}}$ obtained via equation (8).

The development of the dynamic predictive model using a moving window approach combined with Ridge Regression proved to be highly effective. This strategy was designed to capture the evolutionary nature of the biological process and to prevent overfitting that can occur with small data windows. The performance of the final predictive model on the test data was excellent, with an $R^2$ of 0.920 and an MRE of 6.50%. This confirms that the model not only fits the training data well but also generalizes effectively to make accurate short-term predictions on unseen data. The use of Ridge Regression ($\lambda =2$) was crucial in stabilizing the coefficients and ensuring the model's robustness.

Despite the promising results, this study has several limitations that should be acknowledged. First, the model was developed using a dataset from a single anaerobic digestion process spanning 58 days. Its performance and the relative importance of the selected variables may not be directly transferable to other digesters with different feedstock or operating conditions. Second, the initial dataset was expanded through linear interpolation to facilitate the moving window analysis. While this is a valid mathematical procedure, it does not generate new experimental information and could potentially mask high-frequency dynamics not captured by the daily sampling rate.

The analysis of the data from the training and prediction curves is recorded in Table 11. The $R^2$ fit values show a strong correlation (greater than 0.9), and the MRE errors are considerably low (less than 7%), with a comparatively low RMSE for the prediction.

Future research should focus on addressing these limitations. Validating the simplified model on long-term datasets from diverse anaerobic digesters is a critical next step to assess its generalizability. Further investigation into the specific metabolic roles of the top 12 identified microorganisms could provide deeper biological insights to complement the model statistical findings. Finally, exploring non-linear modeling techniques, such as gradient boosting or recurrent neural networks, could potentially capture more complex relationships within the data and further improve predictive accuracy. The ultimate goal would be to integrate such a validated, simplified model with online monitoring sensors for the key microbial populations, paving the way for advanced real-time control and optimization of anaerobic digesters.

4. Conclusions

This study successfully developed a multiple linear regression model to predict methane concentration in anaerobic codigestion using integrated microbial and operational data. The model demonstrated high predictive accuracy (R² = 0.92, MRE = 6.50%) while requiring only 12 key predictors, substantially reducing complexity compared to the initial 31 variable set. Among the relevant findings, the identification of Verrucomicrobia as a significant predictor was particularly noteworthy, as this phylum is typically associated with extreme environments rather than mesophilic digesters, suggesting previously unrecognized ecological adaptations. The overwhelming dominance of microbial indicators over conventional process parameters highlights the critical importance of community dynamics in driving methanogenic performance. Furthermore, the moving window approach with Ridge regularization effectively captured the system's biological evolution while maintaining robustness against overfitting. This modeling approach demonstrates significant potential for practical implementation in rural and resource-limited settings, offering a viable method for methane prediction without sophisticated computational requirements.

Future work should focus on validating this model across diverse reactor configurations and feedstock types to assess its generalizability. Additionally, developing cost-effective molecular monitoring tools for the identified key microbial groups could enable real-time implementation of this predictive approach in practical applications.