Simulation of the biological processes of AD and DF was performed in Aspen Plus V9 software, the Soave-Redlich-Kwong (SRK) method was chosen as it correlates and calculates the mole fractions and activity coefficients. In addition, the SRK model is suitable for vapor-liquid equilibrium (VLE) (Madeira et al., 2017). The starting point for the simulation of AD was the ADM1 model, which relates the stoichiometric equations for the stages of hydrolysis, acidogenesis, acetogenesis, and methanogenesis (Serrano, 2010). The same methodology was also used to simulate the stages of hydrolysis, amino acid degradation, and fermentation for DF, as proposed by Serrano (2011).

For each simulation in Aspen Plus, carbohydrates (e.g. dextrose), complex carbohydrates (e.g. cellulose, hemicellulose), fats (triolein, tripalmitin), proteins, and water mass fractions of the feed substrate were included in the reactor inlet stream. A maximum of 3 U were used in Aspen Plus to simulate all stages of the biological processes (Figure 1).

The first unit is a stoichiometric reactor (RStoic) that simulates the hydrolysis stage for both the AD and DF processes. It calculates the mass flows for carbohydrates, proteins, and fats based on the hydrolysis reactions and fractional conversions (Supplementary Table S1) (Rajendran et al., 2014a).

The second unit is a continuous stirred tank reactor (CSTR) for both the AD and DF processes. The CSTR was selected following the most common operating mode of AD and DF reactors in the literature (Barca et al., 2015; Adarme et al., 2017; Neshat et al., 2017; Rajesh Banu et al., 2020). The CSTR simulates the processes of amino acid degradation, acidogenesis, acetogenesis and methanogenesis for AD, and amino acid degradation and fermentation for DF (Figure 1). A set of 42 reactions was used to describe the main metabolic pathways associated with the bioconversions for the AD and DF processes (Supplementary Table S2) (Rajendran et al., 2014a). The compounds in Table S2 come from the feed substrates and/or from the hydrolysis stage (Supplementary Table S1). All the kinetic reactions are supposed to follow first-order kinetics (Angelidaki et al., 1999). Simulations were performed according to the process simulation model (PSM) proposed by (Rajendran et al., 2014a). PSM model is based on the ADM1, and it includes ammonia inhibition on the calculations of the kinetic constants for the conversion reactions, thus following the model proposed by Angelidaki et al. (1993) and Angelidaki et al. (1999). The model also takes into account variables such as pH and temperature to the calculations of the kinetic constants by correlation functions presented by (Siegrist et al., 1993; Vavilin et al., 1994; Serrano, 2010). All these correlation functions were introduced in the CSTR unit through calculator blocks available in Aspen Plus. Each calculator block had a FORTRAN code to determine the kinetic constants of the conversion reactions. In addition, for DF process a calculator block was integrated into the model to establish the restriction of hydrogen gas production yield from carbohydrates within the values 1–4 mol H₂∙mol hexose⁻¹, where 4 is the maximum stoichiometric yield (Alexandropoulou et al., 2018).

The third unit was developed and applied for the DF process to simulate the homoacetogenic stage, which involves the consumption of hydrogen (H₂) and carbon dioxide (CO₂) gas (Eq. 1) (Saady, 2013). 4 H 2 + 2 C O 2 → C H 3 C O O H + 2 H 2 O                               Δ G ° = 104.0   k J / m o l (1)

H₂ consumption by homoacetogenesis usually increases by increasing the mass transfer limitation of H₂ from the liquid to the gas phase (Saady, 2013). In this study, the minimum conversion rates reported in the literature (Saady, 2013) were considered to simulate low mass transfer limitation conditions in the DF reactors. According to this, homoacetogenesis is carried out in an RStoic block with a conversion factor of 0.1 ± 0.1, as a function of the H₂ and CO₂ composition of the biogas produced during the fermentation step. The methanogenic stage controls the H₂ concentration as H₂ is consumed by hydrogentrophic methanogenesis (pathway n° 39, Supplementary Table S2). Therefore, the homoacetogenic stage is discarded for the AD process.

The residual biomasses selected for this study were coffee mucilage (CFM), cocoa mucilage (CCM) and swine manure (SM), which are related to the most relevant agro-industrial activities in Colombia. Raw samples of CCM, CFM and SM were collected from the departments of Santander and Cundinamarca from the farms on a medium production scale. Then, the total volatile solids (TVS) and total organic load (COD) of the samples were determined by standard methods (APHA, 2005). Total carbohydrates, cellulose, hemicellulose, lignin, and protein were determined by high performance liquid chromatography (HPLC), by using a chromatographer Agilent Model 1200 (Agilent Technology, United States). Total organic carbon (C) was determined by the Walkley Black titration method according to the Colombian technical norm NTC 5167 (2003). Total organic nitrogen (N) and total organic sulfur (S) were determined after acid digestion according to the Kjeldahl method and to the barium sulfate precipitation method, respectively. Then, oxygen (O) and hydrogen (H) contents were estimated stoichiometrically. The results of biomass characterization are summarized in Table 1. These results are in good agreement with those provided by the online database of the Netherlands Energy Research Center (Phyllis2, 2020).

A new category for the substrates had to be created in Aspen Plus due to the database not having this specific kind of complex biomaterials. Therefore, CCM, CFM and SM were defined in the database of Aspen Plus through a set of specific properties, as summarized in Table 1. In the inlet streams for AD and DF models, total fats were divided in equal fractions between triolein and tripalmitin.

The validation of the AD and DF models was carried out by comparison of simulated results to those of batch laboratory tests that aimed at evaluating the biogas production potential of different blends of CFM, CCM, and SM (Cabeza et al., 2016; Hernandez et al., 2018). A Box-Behnken design (Table 2), based on three independent parameters, was used to determine the number of blends for the comparative experiments: 1) the dried mass ratio between CFM and CCM (CFM:CCM ratio), 2) the mass ratio of C to N (C/N ratio), and 3) the total organic load. For each parameter, three levels were considered: the CFM:CCM ratios were 3:1, 1:1, and 1:3; the C/N ratios were 25, 35, and 45; and the organic loads, based on chemical oxygen demand (COD), were 2, 5, and 8 gCOD·l⁻¹. SM was used as N rich substrate to set the C/N ratio of each blend. Then, distilled water was added to set the COD concentrations of each blend.

The experiments were carried out in 250 ml bottles with a total working volume of 200 ml. The reactors were operated under batch conditions and at mesophilic temperature (35°C). The inoculum used for AD and DF experiments was obtained from the biodigester of the wastewater treatment plant of a dairy company (ALPINA, Sopó, Cundinamarca). Before DF tests, a heat shock treatment of the inoculum (30 min at 100°C) was applied to select spore-forming hydrogen-producing bacteria and to destroy hydrogen-consuming archaea that are usually non-spore-forming (Barca et al., 2015). The measurement of methane (CH₄) and hydrogen gas (H₂) was taken daily by volume displacement. NaOH solution (pH > 10) was used as a carbon dioxide gas (CO₂) trap and, therefore, it is assumed that the displaced gas only corresponds to CH₄ in process AD and H₂ in process DF. The duration of the tests was between 10 and 20 days for AD, and between 5 and 12 days for DF experiments. The experimental biogas yields (ml CH₄∙g COD⁻¹ and ml H₂∙g COD⁻¹ for AD and DF, respectively) were calculated based on the cumulated volume of biogas produced from each experiment. All the tests were performed in triplicate and average biogas yields were determined for each experimental condition. As the duration of the tests varied between the different experiments, it was used as the input value of hydraulic retention time for each simulation. All simulations were performed at 35°C and 0.1 MPa. The simulated biogas production yields (ml biogas∙g COD⁻¹) were calculated based on the mass amount of CH₄, H₂, and CO₂, and on the density of each compound at standard conditions (25°C and 0.1 MPa).

The biogas production performances of each substrate blend (Table 2) were further simulated at three different process schemes: anaerobic digestion (AD), dark fermentation (DF), and dark fermentation followed by anaerobic digestion (ADF). The main objective was to identify the most suitable blends and process schemes, for the production of biogas or a specific gas (CH₄, H₂, CO₂), using a comparative simulation approach, thus taking into account relevant factors such as the C/N ratio and the organic load of the feed. Indeed, previous studies have already shown that these factors may significantly affect biogas production performance (Zhang et al., 2014; Xie et al., 2016; Kovalovszki et al., 2017; Mosquera et al., 2020). All the simulations were performed at 35°C and 0.1 MPa, and according to hydraulic retention times of 20 and 10 days for AD and DF, respectively. Then, the biogas yields (ml biogas∙g COD⁻¹) were calculated to assess and compare the performances of the different process scenarios.

The energy recovery potential in Colombia was based on residual biomass production data of 26 departments (out of the 32 Colombian departments) according to the agricultural production reports of the National Department of Statistics, the Colombian Institute of Agriculture, and the National Pork Producers Association (DANE, 2020). A factor of 0.75 was used to estimate the CFM and CCM effectively available (Şenol, 2019), thus assuming 25% of losses in their collection. In addition, the SM availability was calculated as 75% of the total SM from technical farms (ICA, 2018). Once the availability of CFM, CCM, and SM was established, the annual total was divided as daily flow and the ADF process was simulated for each department. The inlet organic load of 26 gCOD·l⁻¹ was used as an input value for all the department simulations based on the COD concentration of the substrates (Table 1) and their total available amounts for each department (Table 3). Temperatures and heat duties, from the results of ADF simulations, were analyzed and gathered to perform the heat integration procedure, according to the pinch method, thus considering the minimum approach temperature ΔTmin of 3°C. Heat transfer from the effluents was integrated into the calculations by considering their cooling down to a temperature of 25°C (Abdou Alio et al., 2021).

For each department, the total daily masses of inlet substrate M_substrate (CCM, CFM, SM) (kg∙day⁻¹), outlet biogas M_biogas (NH₃, CO₂, H₂, CH₄, H₂S) (kg∙day⁻¹), and outlet digestate M_digestate (propionate, butyrate, ethanol, benzene, acetate, furfural, residual carbohydrates, and proteins) (kg∙day⁻¹), were calculated based on the results from simulations. Then, total input power P_input (MJ∙day⁻¹), substrate power P_substrate (MJ∙day⁻¹), biogas power P_biogas (MJ∙day⁻¹), and digestate power P_digestate (MJ∙day⁻¹) were determined by Eqs 2–5, respectively, where LHV is the related lower heating value (MJ/kg) and P_process (MJ∙day⁻¹) represents the energy flow required to heat the feedstock and maintain a constant process temperature in the reactors. P i n p u t = P s u b s t r a t e + P p r o c e s s (2) P s u b s t r a t e = M s u b s t r a t e ⋅ L H V s u b s t r a t e (3) P b i o g a s = M b i o g a s ⋅ L H V b i o g a s (4) P d i g e s t a t e = M d i g e s t a t e ⋅ L H V d i g e s t a t e (5)

The energy recovery potentials E (toe) from biogas and digestate by ADF process were calculated using Eqs 6, 7, respectively, where d represents the days of the year (365) and f_E represents the conversion factor ton of oil equivalent per Joule (2.388ˑ10⁻⁵ toe/MJ). E b i o g a s = P b i o g a s ⋅ d ⋅ f E (6) E d i g e s t a t e = P d i g e s t a t e ⋅ d ⋅ f E (7)

The energy recovery yields for biogas η _biogas (%), digestate η _digestate (%), and total energy recovery η _total (%) were calculated by Eqs 8–10, respectively. η b i o g a s = P b i o g a s P i n p u t (8) η d i g e s t a t e = P d i g e s t a t e P i n p u t (9) η t o t a l = P b i o g a s + P d i g e s t a t e P i n p u t (10)

Finally, a sensitivity analysis of the ADF model was made to assess the effect of nine different organic loads (from 2 to 26 gCOD·l⁻¹) on the total energy recovery yields η _total (%) for the five departments with the highest production of coffee, cocoa and pork (Antioquia, Boyacá, Cundinamarca, Meta, and Santander).

Figure 2 shows the correlation between simulated and experimental biogas production yields for the different blends tested. Overall, the results show an average ratio between simulated and experimental yields (y_sim/y_exp) of 1.05 ± 0.24 for AD and 1.17 ± 0.28 for DF (average values ± standard deviations), respectively, thus indicating a good general agreement between the data.

With the only exceptions of blends 1, 2, and 13, the ratio y_sim/y_exp for AD varied from 0.79 to 1.26 (Figure 2A), thus indicating that the AD model gives a good estimation of CH₄ production yields. The lowest differences were observed for blends 3, 7, and 10 (y_sim/y_exp ratios of 0.92, 1.00, and 1.09, respectively), whereas the highest differences were observed for blends 2 and 13 (y_sim/y_exp ratios of 0.70 and 1.79, respectively). The results indicate that the differences between simulated and experimental data for AD increase by decreasing the organic loads and/or by increasing the SM content of the blends. Higher organic loads probably resulted in dominant organic compounds in the blends used for the experiments, thus leading to more stable and more predictable conversion pathways. Moreover, variability in the single substrate composition used for the different experiments, especially for SM, may further explain the differences between simulated and experimental yields.

With the only exceptions of blends 4, 7, and 12, the ratio y_sim/y_exp for DF varied from 0.77 to 1.28 (Figure 2B), thus confirming the validity of the DF model. The lowest differences were observed for blends 3 and 6 (y_sim/y_exp ratios of 1.03 and 1.06, respectively), whereas the highest differences were observed for blends 4, 7, and 12 (y_sim/y_exp ratios of 1.72, 1,94, and 1.43, respectively), thus indicating an overestimation of H₂ production yields compared to the experimental data. The results of DF seem to indicate that the differences between simulated and experimental yields increase by increasing the CCM and/or CFM loads of the blends. Most probably, incomplete hydrolysis of complex carbohydrates (e.g. cellulose, hemicellulose) may account for the lower experimental yields for blends at high loads of CCM and/or CFM (blends 4 and 7, respectively).

The results of the validation model suggest that variability on single substrate composition between the different experiments and incomplete hydrolysis of carbohydrates may account for the differences between simulated and experimental yields. However, using the same AD and DF models for evaluating the potential of biogas production potential from model substrates may overcome these technical limitations of laboratory experiments and thus allow comparison of the potential of energy recovery potentials from different blends and different process schemes.

The total biogas production yields and compositions obtained by the simulations of AD, DF, and DF followed by AD (ADF) processes, for the different substrate blends (Table 2), are summarized in Figure 3.

The biogas production yields for DF varied between 147 and 405 ml biogas·gCOD⁻¹, and the H₂ and CO₂ content of biogas varied between 65–70% and 30–35%, respectively (Figure 3C). These results are in good agreement with the theoretical compositions of biogas from the fermentative pathways used for DF simulation (n° 39 and 40, Supplementary Table S2) (Gomes et al., 2015). The lowest biogas production yield for DF was observed for blend 9, where the dominant substrate is SM (more than 80% based on total organic load) and COD concentration is 2 gCOD∙l⁻¹. Instead, the highest biogas production yield for DF was observed for blend 7, where the dominant substrate is CFM and COD concentration is 5 gCOD∙l⁻¹. These differences in biogas production yields are probably primarily related to the different compositions of proteins, fats, and carbohydrates considered for CFM and SM (Table 1). Higher carbohydrate content of CFM resulted in higher H₂ production yields for blend 7, as carbohydrates represent the most suitable substrate for DF (pathways n° 39 and 40, Supplementary Table S2), whereas higher protein content of SM resulted to lower C/N ratios that may inhibit DF as described by the ammonia inhibition modules integrated into the model (Madeira et al., 2017; Serrano, 2011). These results are in good agreement with the findings of previous experimental studies (Zhu et al., 2008; Ghimire et al., 2015), which have indicated that low C/N ratio and high concentration of ammonium ion (NH₄ ⁺) inhibit H₂ production during DF when SM is used as the only substrate, with maximum production yields of 3.65 ml H₂∙VS⁻¹. Indeed, higher production yields (from 68 to 652 ml H₂∙gCOD⁻¹) were obtained by co-fermentation of substrate blends (e.g. swine manure, sewage sludge, food waste, agro-industrial waste) (Kim et al., 2004; Perera et al., 2012; Hernández et al., 2014).

The biogas production yields for AD varied between 295 and 459 ml biogas·gCOD⁻¹, thus showing a lower variability according to the different blends rather than DF. The CH₄, CO₂ and H₂ biogas components were 57–70%, 29–42%, and 0.3–0.4%, respectively (Figure 3B). Higher yields of biogas production were observed for blends containing higher amounts of CFM and CCM rather than SM. Overall, biogas production yields obtained by simulations are within the ranges reported in the literature (61–650 ml CH₄∙gCOD⁻¹). This large variability in experimental yields is mainly attributed to the different organic substrates and operating parameters (e.g. retention time, temperature, inoculum) between the different studies (Garfí et al., 2011; Hernández and Rodríguez, 2013; Astals et al., 2015). The blends 12 and 13, with equal CFM and CCM contents, showed the highest biogas yields (433 and 459 ml biogas·gCOD⁻¹, respectively) with the highest CH₄ content (70 and 69%, respectively). Therefore, they appear to be the most suitable blends for AD. These results are consistent with the experimental results of Wang et al. (Wang et al., 2018), who indicated that biogas production performances of AD are significantly affected by the ratios of swine manure, corn stove, and cucumber residues in the co-digested mixture.

The biogas production yields for DF followed by AD (ADF) varied between 602 and 864 ml biogas·gCOD⁻¹, and the CH₄, CO₂, and H₂ content of biogas varied between 25–51%, 30–35%, and 17–40%, respectively (Figure 3D). The highest total biogas production yields were observed for blends 7 and 8 (863 and 864 ml biogas·gCOD⁻¹, respectively), where the CFM and CCM are above 65% of the total organic load. Additionally, one of the substrates, CFM or CCM, dominates over the others. The lowest yields (<160 ml biogas·gCOD⁻¹) were obtained for blends 5, 6, 9, and 10, where SM is greater than 80% over the total organic load. However, despite the lowest biogas production yields, these blends also show the highest CH₄ contents in biogas. This appears to confirm that SM is more suitable for CH₄ than for H₂ production, most probably because it has a lower content in carbohydrates compared to the other substrates. While, blends with an SM load equal or lower than those of CFM and CCM (blends 7, 8, 11, 12, and 13) generate biogas with CH₄ and H₂ contents comprised between 28 and 40%. ADF showed higher yields of total biogas production compared to AD and DF, this for all the substrate blends studied. As shown in Figure 3A, total biogas production yields for ADF were 1.3–2.3 times higher than DF, and 1.5–2.4 times higher than AD, thus indicating a greater conversion efficiency of organic substrate to biogas. This is coherent with the assumptions of the ADM1 model that was used for the simulations, as organic products from DF (e.g. volatile fatty acids) represent viable substrates for methanogenesis in the AD process (Hernández et al., 2014), thus increasing the total biogas production yields (Yeshanew et al., 2018).

The results of the simulations confirm the effectiveness of co-digestion of CFM, CCM, and SM, in a two-step process of DF followed by AD compared to the individual process steps. The higher content in proteins, fats, and hemicellulose of SM supports the production of fatty acids, amino acids, and carbohydrates by hydrolysis in the DF step (Supplementary Table S1). These byproducts are further converted and used as substrates for CH₄ production by acidogenesis, acetogenesis, and methanogenesis in the AD step (Supplementary Table S2). Instead, the higher content in carbohydrates of CFM and CCM conducted to the production of acetate and H₂ by fermentation in the DF step (Supplementary Table S2), which are used as substrates for CH₄ production by acetotrophic and hydrogenotrophic methanogenesis in the AD step (pathways 37 and 38, Supplementary Table S2). Moreover, the simulation results show that different organic loadings do not seem to have a significant impact on the production yields and biogas composition for the different blends and processes studied. This behavior was related to all blends and processes simulated with the same stoichiometric ratios, considering residence times long enough to reach conversion equilibrium and without considering excess substrate inhibition. However, despite showing similar production yields, higher organic loads may give higher flow rates of biogas production in real operating conditions, as more organic substrate is available for biogas production per unit of time (Barca et al., 2015).

Among the 26 departments considered in this study, 17 report coffee production more than twice higher than cocoa production (DANE, 2020). Only eight departments (Bolívar, Córdoba, La Guajira, Meta, Arauca, Putumayo, and Vichada) show higher cocoa production than coffee production, and their total coffee production amounts are relatively low (<210 tons of coffee/year) compared to the other departments (>2000 tons of coffee/year). According to DANE data (2017), the departments of Huila, Cauca, Tolima, Caldas, and Antioquia generates 73% of the Colombian coffee production and they are therefore considered as the main Colombian coffee region. Additionally, their coffee production is more than 15 times higher than their cocoa production, thus indicating a much lower availability of CCM rather than CFM (Supplementary Figures S2, S3).

Table 3 summarizes the CFM, CCM, and SM availability data for 26 Colombian departments in 2017. Overall, the data indicates that CFM production dominates over CCM production in the center and south of the Andean region (Antioquia, Huila, Tolima, Caldas, Risaralda, Quindío), in the Pacific region (Chocó, Valle del Cauca, Nariño, Cauca), and more generally in the south-west of Colombia, whereas CCM production dominates over CFM production in the Orinoquía region (Arauca, Casanare, Meta, Vichada), in the Caribbean region (Atlántico, Bolívar, Córdoba, Sucre, La Guajira Cesar, Magdalena), in the north of the Andean region (Boyacá, Norte de Santander, Cundinamarca, Santander), and more generally in the north-east of Colombia. The Amazon region does not show significant productions of CFM and CCM, with the only exception of the department of Caquetá. SM availability per department in 2017 varies from 3,347 ton (Vichada) to 1084000 ton (Antioquia). The mass ratio of available SM over the sum of available CFM, CCM, and SM (SM/total) varies from 0.878 (Santander) to 1.000 (Atlántico), thus indicating a larger availability of SM compared to CFM and CCM for all the departments.

For each department, the potential recovery of biogas and digestate was evaluated by comparative ADF simulations according to the total available amounts of CFM, CCM, and SM in 2017, and considering for all the simulations an inlet COD concentration of 26 g COD∙l⁻¹. As shown in Table 4, the total mass of biogas recovered per department varied from 307.6 ton (Vichada) to 102200 ton (Antioquia), whereas the total mass of digestate recovered per department varied from 146.3 ton (Vichada) to 47,580 ton (Antioquia). The main components of biogas were CO₂ and CH₄, with an evaluated production potential ranging between 176.6–58,960 ton and 105.7–34,950 ton, respectively, whereas the main components of digestate were ethyl cianoacetate and propionate, with an evaluated production potential ranging between 11.19–3,693 ton and 10.08–3,295 ton, respectively.

Table 5 summarizes the energy recovery potential (toe) from the biogas and digestate obtained for each department based on production data of 2017, as well as the energy recovery yields for biogas η _biogas (%), digestate η _digestate (%), and the total energy recovery yields η _total (%) (Equations 8–10). The highest energy recovery potentials from biogas were observed for the departments of Antioquia, Córdoba, and Cundinamarca (48,990, 10,160, and 13,460 toe, respectively), which are the departments with the highest cumulated productions of CFM, CCM, and SM. Instead, the lowest energy recovery potentials from biogas were observed for the departments with the lowest cumulated productions of CFM, CCM, and SM (Putumayo and Vichada, 866 and 148 toe, respectively). According to the energy consumption data in 2017 (IEA, 2020), the total energy recovery potential from the biogas (151.6 ktoe) could supply up to 70% of the energy demand for agriculture/forestry activities in Colombia (215 ktoe for 2017).

As shown in Table 5, the energy recovery yields for biogas η _biogas (%) varied between 34.3 and 39.3%. The highest η _biogas (%) were observed for the departments of Cauca, Huila, Santander, and Tolima (36.3, 36.3, 39.3, and 36.5%, respectively), which are the departments with the highest CFM/SM and/or CCM/SM production ratios. On the contrary, the departments with lower CFM/SM and/or CCM/SM production ratios (Atlántico, Córdoba, Sucre, and Vichada) usually present lower energy recovery yields for biogas η _biogas (%) and higher energy recovery yields for digestate η _digestate (%). This seems to confirm that although SM is less effective than CFM and CCM in biogas production, it is more effective for recovering other by-products from digestate (e.g. propionate, butyrate, ethanol, benzene, acetate, and furfural). As shown in Table 5, the total energy recovery yields η _total (%), which include biogas and digestate energy (Equation 10), varied between 48.6% (Atlántico, Sucre, and Vichada) to 52.5% (Santander).

The results show higher total energy recovery yields η _total (%) for the departments with higher CFM/SM and/or CCM/SM production ratios, confirming that CFM and CCM are more effective than SM for energy recovery by ADF processes. It should be noted, however, that since this is a theoretical study of energy recovery potential, factors such as seasonal variations in agro-industrial production and variability in substrate composition between different production sites due to breeding and cropping techniques, as well as soil and climate characteristics, are not considered.

One of the impact factors on biogas production yields is the input load (Jiang et al., 2020). Figure 4 shows the evolution of total energy recovery yield η _total (%) as a function of inlet organic load for the five departments with the highest production of CFM, CCM, and SM (Antioquia, Boyacá, Cundinamarca, Meta, and Santander). The η _total (%) values in Figure 4 were obtained by comparative ADF simulations according to the CFM, CCM, and SM production data in 2017. The results indicate that η _total (%) increases from 10.2–11.4% to 48.4–51.4% by increasing inlet COD from 2 to 26 gCOD·l⁻¹ for all the departments. This was primarily related to the lower energy consumption for the initial heating of the feed, since the higher the COD concentration, the lower the feed flow to be heated. However, the rate of η _total (%) increase appears to decrease by increasing organic load, thus suggesting an optimal inlet COD concentration to obtain higher η _total (%) while maintaining acceptable yields of substrate conversion. Moreover, increasing the inlet COD concentration over the average value of the available raw mixtures (26 gCOD·l⁻¹) would require further energy consumption (e.g. drying and membrane processes) that may affect the energy recovery efficiency of the full process. Among the departments in Figure 4, the department of Santander shows the highest η _total (%) for all the inlet COD concentrations. This was probably because Santander has the highest ratios (CFM + CCM)/SM and C/N, thus resulting in higher biogas production performances compared to the other departments.

In recent years the Colombian economy has come to strongly depend on hydrocarbon production from non-renewable natural resources. In April 2020, 16.6 million barrels of crude oil were exported, which represented a 19.7% drop compared to April 2019 (DANE, 2020). Biofuel production from alternative renewable sources such as agro-industrial wastes may represent a sustainable way to overcome the shortage of natural deposits of fossil carbon. During the last 5 yr, the increase in pork breeding, coffee, and cocoa production (DANE, 2020) suggests an interesting potential market of SM, CFM, and CCM valorization as organic substrates for biogas production. However, the fluctuation in production of these residual biomasses may represent a weakness for perspective evaluations, and therefore a regular updating of inlet data for the models is recommended to better follow future production dynamics. It should also be noticed that most of the coffee, cocoa, and pork production activities are located in rural areas, and, currently, more than 20% of transport roads in Colombia are without pavement (DANE, 2020). This may significantly hinder the collection of co-substrates from different production areas, and it may also generate supplementary transport costs. The bioprocess simulation approach developed in this study may be followed to evaluate energy recovery potential from other types of locally available substrates (IEA, 2020), thus showing flexibility of adaptation to local contexts. Furthermore, the integration of thermochemical (e.g. gasification, pyrolysis, hydrothermal liquefaction) and separation (e.g. membrane, extraction, distillation) processes may allow the recovery of further by-products (e.g. nutrients, biochar, amino acids, organic acids) from the effluents of AD and DF processes (Cebreiros et al., 2017; Arora et al., 2018; Rajesh Banu et al., 2020; Uddin et al., 2021). Therefore, an accurate life cycle analysis (LCA), integrating energy balances and operating costs for each unit of the process, would be required to identify the best design, location, and size of the biorefinery plants in Colombia.