<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamic modeling of biological methanation for different reactor configurations: An extension of the anaerobic digestion model No. 1

Juan Acosta-Pavas, Carlos Robles-Rodríguez, Jérôme Morchain, Claire Dumas, Arnaud Cockx, César Aceves-Lara

► To cite this version:

Juan Acosta-Pavas, Carlos Robles-Rodríguez, Jérôme Morchain, Claire Dumas, Arnaud Cockx, et al. Dynamic modeling of biological methanation for different reactor configurations: An extension of the anaerobic digestion model No. 1. Fuel, 2023, 344, pp.128106. 10.1016/j.fuel.2023.128106. hal-04184348

HAL Id: hal-04184348

https://hal.science/hal-04184348v1

Submitted on 31 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
1 Dynamic Modeling of Biological Methanation for Different Reactor

2 Configurations: An Extension of the Anaerobic Digestion Model No. 1

3 Juan C. Acosta-Pavasª*, Carlos. E Robles-Rodríguezª, Jérôme Morchainª, Claire Dumasª, Arnaud

4 Cockxª, César A. Aceves-Laraª

5

6 ªTBI, Université de Toulouse, CNRS, INRAE, INSA, Toulouse, France

7

8 Corresponding author: acostapa@insa-toulouse.fr (Juan C. Acosta-Pavas)

9

10 Abstract

11 In biological methanation, the methane produced by anaerobic digestion (AD) is upgraded with the
12 addition of syngas. Several mathematical models have been developed to represent the AD process.

13 However, the modeling of biomethanation is still under development. In this work, an extension of the
14 anaerobic digestion model (ADM1_ME) was proposed to describe the dynamics of biomethanation.

15 The model considered the addition of syngas flow rate (H₂/CO) and the adaptation of the volumetric
16 mass transfer coefficient for two different reactor configurations: bubble column reactor and
17 continuous stirred tank reactor operating at mesophilic and thermophilic conditions. A sensitivity
18 analysis using the Sobol' and Morris methods was performed for this model, where 14 parameters
19 were selected for model calibration. Simulation results showed an accurate fit for two experimental
20 operating conditions from literature with a RMSE <5.15. The results showed the feasibility of the
21 ADM1_ME to describe the biomethanation process at different operational conditions and reactor
22 configurations.

23

24 **Keywords:** Anaerobic digestion model, Biological methanation, Sensitivity analysis, Parameter
25 estimation

26
27

# 1. Introduction

28

29

30

31

The production of biogas (a mixture of CH₄ and CO₂) by anaerobic digestion (AD) is currently one of the most promising options in terms of bioenergy production [1]. This biogas can be used locally to generate heat or electricity without any additional processing. However, it needs to be refined or purified to be used as vehicle fuels or for injection into the gas grid system [2,3].

32

33

34

35

36

37

38

Biogas could endure post-treatment alternatives either to remove impurities (e.g., H₂S, excess of water) or to be upgraded into biomethane (95 – 99 % CH₄) for a further injection in the gas grid.

Biogas upgrading consists of increasing the CH₄ concentration by removing CO₂ [4,5]. Some of the technologies employed in biogas upgrading include physical absorption (water, amine, and organic scrubbing), pressure swing adsorption, and membrane separation [2,4]. Another process that has gained interest and is currently under development is biomethanation, also called biological methanation [6].

39

40

In the biological methanation process the CO₂ contained in the biogas is converted into CH₄ by using hydrogen coming from the addition of syngas, which is generally composed of CO₂, CO, and H₂ [5].

41

42

43

44

45

Biomethanation involves a complex microbial consortium, whose composition changes due to operational conditions, such as temperature, pH, hydraulic retention time, and syngas composition [7,8]. Biomethanation can be performed either *in-situ* where syngas is introduced into the AD reactor, or *ex-situ* where a microbial consortium, coming from the AD and adapted for H₂ consumption, is introduced into a second bioreactor to convert syngas into a high-purity CH₄ [5].

46

47

Until now, the development of dynamic models for AD has been focused on the incorporation of detailed knowledge of the process stages, the microorganisms involved, and the operating conditions.

48

49

50

51

52

53

One of the first dynamic models to investigate AD was proposed by Mosey [9], who developed a model to consider how the microorganisms managed to control the pH value and the redox potential of their growth medium. A dynamic model describing AD from several types of wastes was developed in [10]. In their work, the substrate was described by its composition in terms of basic organic components, i.e., carbohydrates, lipids, proteins, volatile fatty acids, long-chain fatty acids, and inorganic components. This partition allowed the authors to simulate the dynamic changes during the
process with different types of substrates. Batstone et al. [11] proposed the well-known Anaerobic
Digestion Model No. 1 (ADM1), which included multiple steps describing the biochemical and
physicochemical reactions involved in AD. Regarding the biochemical reaction, the ADM1 considers
the disintegration of components, such as carbohydrates, lipids, and proteins into particulate
constituents that are further hydrolyzed into soluble monomers which are precursors of CH₄ formation.
The physicochemical reactions describe ion associations and dissociations, as well as the gas/liquid
transfer phenomena. The ADM1 has been modified [12] to solve stiffness problems, especially due to
the mass transfer equations. Thus, these authors proposed a Benchmark Simulation Model (BSM).
Several researchers have adapted the BSM to consider inhibition by free ammonia in high-solid sludge
fermentation [13], to design optimal continuous operation of experimental anaerobic digestion [14], or
to simulate the dynamic behavior of a pilot-scale process for two-stage anaerobic digestion of sewage
sludge [15].

Modeling of biomethanation has been scarcely studied. To our knowledge, only two works have been
reported. Grimalt-Alemany et al. [7] proposed two structured models to describe the mesophilic and
the thermophilic syngas biomethanation processes in batch mode. Each model presented a different
structure based on catabolic routes as a function of the operating conditions. All biomass growth
processes were made thermodynamically consistent by including a thermodynamic potential factor.
Although these adaptations improved the predictive capacity of the models, the studied carbon source
was only limited to the added gas which could hamper any straightforward adaptation to other sources,
such as agro-industrial residues, sludge, or sugars. Additionally, the different volumetric mass transfer
coefficients were calculated experimentally under specific mesophilic and thermophilic conditions.
Although good, the validation of results and the calculation of the mass transfer coefficient from
experiments make the adaptation of the model to other reactor configurations and conditions difficult.
Sun et al. [16] proposed an extension of the ADM1 considering biochemical reactions for the CO
contained in the syngas under mesophilic conditions in continuous mode. In their work, a volumetric
mass transfer equation was developed from the two-film theory to describe the mass transfer process.
This consideration makes it difficult to apply the model to different process conditions and increases
the number of parameters in the process.

Both models are a good basis for modeling the biomethanation process. However, they present
deficiencies in considering the type of reactor, which directly affects how the mass transfer
phenomenon is represented. This is highly important since H₂ transfer is one of the limiting factors of
biomethanation [17]. Additionally, the model from Sun et al. [16] does not present the dynamic
behavior of components, such as CO₂, which complicates the comprehensive analysis and the closing
of the carbon balance in the biogas produced from AD.

This work aims to generate advances in the understanding of the dynamics of the biomethanation
process by extending the Anaerobic Digestion Model No. 1 (ADM1_ME) to consider in-situ syngas
addition. The model could be adapted to different substrates, agro-industrial waste, sludge, and sugars,
and considers the addition of syngas at mesophilic and thermophilic conditions. The ADM1_ME
represents the dynamic behavior of CH₄, H₂, CO, and CO₂ in liquid and gas phases, the inclusion of
CO as a substrate of the process, and the adaptation of the volumetric mass transfer coefficient for two
different configurations of bioreactor: bubble column reactor (BCR) and continuous stirred tank
reactor (CSTR). This new model is based on the need to have a generic model for the biomethanation
process with the capacity to be adapted to mesophilic and thermophilic conditions, to different
substrates, and to allow the dynamic analysis of the volumetric mass transfer coefficients for different
types of reactors.

## 2. Description of the Biomethanation Process

Biomethanation is a process in which the biogas produced through the well-known AD is upgraded by
the biological conversion of CO₂ and syngas to obtain high-purity CH₄ [6]. In AD the organic matter,
such as agricultural residues, organic effluents from the food industry, animal manure, or
waste/wastewater residues are transformed through the synergistic work of a variety of
microorganisms into a mixture of CH₄ and CO₂ through four steps: (i) hydrolysis, (ii) acidogenesis,
(iii) acetogenesis and (iv) methanogenesis [18]. The biogas produced in the AD contains between 50 -
75% of CH₄, 25 - 50% of CO₂, and 2-7% water vapor [3,4,19]. Hereby, the process is extended to
biomethanation which includes CH₄ production from a gas load converting an inlet flow of H₂ and CO
into high-purity CH₄ [16].
107
108

**Hydrolysis.** In this step, the fermentative bacteria excrete enzymes that transform complex organic
polymers (carbohydrates, proteins, and lipids) into soluble monomers, such as monosaccharides,
amino acids, and long-chain fatty acids (LCFA). This process generally takes place on the surface or
near acidogenic bacteria as it involves exo-enzymes secreted by hydrolytic bacteria, such as
Clostridia, Bacteroides, Fusobacterium, Butyrivibrio, Micrococci, Streptococcus, and Selenomonas
[20,21].
109
110
111
112
113
114

**Acidogenesis.** Throughout the acidogenesis, the dissolved monomers or oligomers, amino acids, LCFA, in general, the components produced in the hydrolysis step undergo degradation reaction. These components are diffused into the acidogenic bacteria through the cell membrane and later fermented or anaerobically oxidized [22] and produce mainly volatile fatty acids (VFA), such as propionate, butyrate, valerate, acetate, as well as new cell material. This step is carried out by the action of bacteria of the genera *Bacillus* sp., *Pseudomonas* sp., *Clostridium* sp., and *Bifidobacterium* sp. [18,20].
115
116
117
118
119
120
121

**Acetogenesis.** The VFA produced in the acidogenesis are reduced and transformed into acetate, H₂, and CO₂, as well as in new cellular material by the action of bacteria of the genera *Clostridium*, *Syntrophomonas* sp., *Syntrophobacter* sp. [20,21].
122
123
124

**Methanogenesis.** Acetate and H₂, are converted into CH₄, and CO₂, as well as in new cellular material
in two types of processes, hydrogenotrophic methanogenesis and acetoclastic methanogenesis by the
strictly anaerobic methanogens of the phyla Euryarchaeota: *Methanobacteriales*, *Methanococcales*,
*Methanomicrobiales*, *Methanosarcinales*, and *Methanocellales*.
125
126
127
128

In hydrogenotrophic methanogenesis, the CO₂ is reduced into CH₄ using H₂ as a reduction agent [23]
by the action of hydrogenotrophic methanogens, *Methanobacterium*, *Methanogenium*,
*Methanocorpusculum*, *Methanobrevibacter*, *Methanoculleus*, *Methanospirillum*,
*Methanothermobacter*, *Methanococcoides*, *Methanomethylovorans*, *Methanosarcina*. In acetoclastic
methanogenesis, the acetate is decarboxylated and converted into CH₄ by the action of acetoclastic
129
130
131
132
133
methanogens, *Methanosaeta*, *Methanococcoides*, *Methanomethylovorans*, and *Methanosarcina* [18,20,22,24–26].

**Biomethanation.** Biomethanation occurs in the last stage of AD. The syngas loading (commonly a combination of H₂/CO/CO₂) can be used to improve the process and convert the H₂ and CO₂ into CH₄ [5]. The H₂ is transformed by the hydrogenotrophic methanogens with CO₂ consumption. Although this is a well-known route, CO consumption is still unclear [16]. The CO can be transformed indirectly into H₂ by carboxydotrophic hydrogenogenesis (water gas shift), then into acetate by CO-acetogenesis, or CO-homoacetogenesis, and finally transformed into CH₄ through the hydrogenotrophic and acetoclastic methanogenesis [27]. Table 1 presents a summary of the reactions involved in the biomethanation process.

<table><caption>**Table 1.** Reactions involved in biological methanation. Adapted from [6,16,23,28–30].</caption><thead><tr><th>Reactions</th><th>ΔG°'(kJ/mol)</th></tr></thead><tbody><tr><td><strong>Acidogenesis reactions:</strong></td><td></td></tr><tr><td>Acetate: C<sub>6</sub>H<sub>12</sub>O<sub>6</sub> + 2H<sub>2</sub>O → 2CH<sub>3</sub>COOH + 4H<sub>2</sub> + 2CO<sub>2</sub></td><td>-206</td></tr><tr><td>Butyrate: C<sub>6</sub>H<sub>12</sub>O<sub>6</sub> → 2CH<sub>3</sub>CH<sub>2</sub>CH<sub>2</sub>COOH + 2H<sub>2</sub> + 2CO<sub>2</sub></td><td>-254</td></tr><tr><td>Propionate: C<sub>6</sub>H<sub>12</sub>O<sub>6</sub> + 2H<sub>2</sub> → 2CH<sub>3</sub>CH<sub>2</sub>COOH + 2H<sub>2</sub>O</td><td>-279.4</td></tr><tr><td><strong>Acetogenesis reactions:</strong></td><td></td></tr><tr><td>Propionate: CH<sub>3</sub>CH<sub>2</sub>COOH + 2H<sub>2</sub>O → CH<sub>3</sub>COOH + 3H<sub>2</sub> + CO<sub>2</sub></td><td>+76</td></tr><tr><td>Butyrate: CH<sub>3</sub>CH<sub>2</sub>CH<sub>2</sub>COOH + 2H<sub>2</sub>O → 2CH<sub>3</sub>COOH + 2H<sub>2</sub></td><td>+48.4</td></tr><tr><td><strong>Methanogenesis reactions:</strong></td><td></td></tr><tr><td>Hydrogen: 4H<sub>2</sub> + CO<sub>2</sub> → CH<sub>4</sub> + 2H<sub>2</sub>O</td><td>-130.7</td></tr><tr><td>Acetate: CH<sub>3</sub>COOH → CH<sub>4</sub> + CO<sub>2</sub></td><td>-31.0</td></tr><tr><td><strong>Carbon monoxide reactions:</strong></td><td></td></tr><tr><td>4CO + 2H<sub>2</sub>O → CH<sub>3</sub>COOH + 2CO<sub>2</sub></td><td>-176</td></tr><tr><td>2CO + 2H<sub>2</sub> → CH<sub>3</sub>COOH</td><td>-67</td></tr><tr><td>CO + H<sub>2</sub>O → CO<sub>2</sub> + H<sub>2</sub></td><td>-20</td></tr></tbody></table>

The biomethanation process can be performed in two ways: *in-situ* (directly in the AD reactor), or *ex-situ* (in a separate unit). The advantage of *in-situ* biomethanation is the reduction of infrastructure costs, due to the use of only one reactor where syngas is directly added allowing AD and biomethanation to occur simultaneously. The main drawback of *in-situ* biomethanation is that AD can be affected by the high concentration of gas. For instance, the H₂ added into the system can inhibit the early stages of the process. On the other hand, *ex-situ* biomethanation takes place in a separate external
reactor, typically adapted to suit the hydrogenotrophic methanogens. Regardless of the type of configuration, microorganisms, and metabolic pathways are similar [5,6,31].

**3. Anaerobic Digestion Model Extension (ADM1_ME)**

In this work an extension of the Anaerobic Digestion Model No. 1 (ADM1_ME) is proposed based on the biomethanation scheme presented in Fig. 1. The scheme includes acidogenesis, acetogenesis, methanogenesis and its extension to consider CO and $H_2$ addition.

<figure><img src="image_2.png" /></figure>

**Fig. 1.** Biomethanation scheme implemented in the ADM1_ME including: ($\mu_{su}$) acidogenesis from sugars, ($\mu_{pro}$) acetogenesis from propionate, ($\mu_{bu}$) acetogenesis from butyrate, ($\mu_{ac}$) acetoclastic methanogenesis, ($\mu_{H_2}$) hydrogenotrophic methanogenesis, and ($\mu_{CO}$) acetogenesis and carboxydotrophic hydrogenogenesis from carbon monoxide. $I_{H_2,k}$ concerns inhibition with respect to $H_2$, and $I_{CO,k}$ concerns inhibition with respect to CO.

**3.1. Mass Balances**

Mass balances are derived from the ADM1 proposed by Batstone et al. [11], which considers disintegration, hydrolysis, and uptake of the various components as well as biomass decay. Differently from ADM1, the model extension ADM1_ME considers the uptake of sugar, volatile fatty acids (e.g., butyrate, propionate, and acetate), the uptake of $H_2$ and CO, and the decay of biomass. The model
describes three types of variables: soluble ($S_{liq,j}$), particulate ($X_k$) and gas ($S_{gas,i}$) components. Particulate components are considered to be part of the biomass as in the ADM1 (see supplementary material Eqs. (S.1) - (S.34) for more details). The model is here rewritten as Eqs. (1) – (3).

**State variables in the liquid phase ($S_{liq,j}$):**

$$ \frac{dS_{liq,j}}{dt} = \frac{q_{liq}^{in}}{V_{liq}}(S_{liq,j}^{in}-S_{liq,j}) + \sum_{k} Y_{k}f_{j,k}\mu_{k} - N_{i} \quad (1) $$

**State variables in biomass ($X_k$):**

$$ \frac{dX_k}{dt} = \frac{q_{liq}^{in}}{V_{liq}}(X_k^{in}-X_k) + Y_k\mu_k - \mu_{k,dec} \quad (2) $$

**State variables in the gas phase ($S_{gas,i}$):**

$$ \frac{dS_{gas,i}}{dt} = \frac{q_{gas}^{in}}{V_{gas}}S_{gas,i}^{in} + N_i \left( \frac{V_{liq}}{V_{gas}} \right) - \frac{q_{gas}}{V_{gas}}S_{gas,i} \quad (3) $$

Sub-index $j \in [1,8]$ denotes glucose, butyrate, propionate, acetate, $H_2$, $CH_4$, CO, and $CO_2$ in the liquid phase. The $H_2$, $CH_4$, and CO are expressed in $gCOD/L$, and $CO_2$ is expressed in $mol/L$. Chemical Oxygen Demand (COD) is the amount of oxygen needed to degrade the organic matter into $CO_2$ and $H_2O$. It is important to mention that $CO_2$ would be expressed in *moles* instead of *COD*, as suggested in [11]. Sub-index $k \in [1,6]$ reads for the biomass that degrade glucose, butyrate, propionate, acetate, $H_2$, and CO, respectively. For the gas phase, the sub-index $i \in [1,4]$ corresponds to $H_2$, $CH_4$, CO and $CO_2$. The inlet flow rates of liquid and gas are represented by $q_{liq}^{in}$ and $q_{gas}^{in}$, respectively, while $q_{gas}$ denotes the outlet gas flow rate. $V_{liq}$ and $V_{gas}$ are the liquid and gas volumes, respectively, $S_{liq,j}^{in}$, $S_{gas,i}^{in}$, and $X_k^{in}$ hold for the inlet concentration of the component $j$ in the liquid phase, the inlet concentration of component $i$ in gas phase, and the inlet concentration of biomass $k$ in the liquid phase. $Y_k$ is the yield of biomass $k$, $f_{j,k}$ the stoichiometric coefficients; $\mu_k$ and $\mu_{k,dec}$ the growth and decay rate of biomass $k$, and $N_i$ the mass transfer rate of component $i$.

In Eq. (1), the term $-\sum_{k=1}^{6}(\sum_{j=1}^{7} C_j v_{j,k} \mu_k)$ is introduced in the $CO_2$ mass balance for the liquid phase (Eq. (S.8) supplementary material), based on the original version of the ADM1 [11] (without considering the components associated with the hydrolysis). This term consists of describing the
fractionation of inorganic carbon, the composition of various species as well as a standard biomass composition [12].

## 3.2. Mass Transfer Rate Definition

To relate the liquid and gas balances, the gas-liquid mass transfer rates are expressed as in Eq. (4).

$$N_i = k_L a_i (S_{liq,j} - \gamma_{COD,i} H_i P_{gas,i}) \quad (4)$$

where $N_i$ is the flux of species H₂, CH₄, and CO, expressed as gCOD/L/d, and CO₂, expressed as mol/L/d. $k_L a_i$ is the volumetric mass transfer coefficient of component i, and $(S_{liq,j} - \gamma_{COD,i} H_i P_{gas,i})$ is the driving force. $H_i$, and $P_{gas,i}$ are Henry's law equilibrium constant, and partial pressure of component i, respectively. $\gamma_{COD,i}$ is a conversion factor between the moles and gCOD of a component i, e.g., 16 gCOD/mol for H₂ and CO, and 64 gCOD/mol for CH₄.

The volumetric mass transfer coefficient ($k_L a$) was set constant (200 1/d) in the ADM1 [11]. However, for the biomethanation process, it is necessary to integrate the effect that the addition of gas has on mass transfer. Depending on the reactor configuration, well-established correlations for BCR Eq. (5) and CSTR Eq. (6) have been used to estimate $k_L a$. Both rely on the superficial gas velocity defined as the inlet gas flow rate divided by the cross-sectional area of the vessel ($U_G = q_{gas}^{in}/A$). If there is no addition of gas to the system, $U_G$ is zero. In a BCR, the sparger type is crucial to ensure efficient gas distribution [32]. For a CSTR, the energy injected in the system controls the bubble size and turbulence, and thus correlates with the mass transfer coefficient. This energy dissipation per unit of volume in the presence of gas, the so-called gassed power input (Power), depends on the impeller type, rotation speed and aeration number [33,34].

**Bubble column reactor (BCR).** The volumetric mass transfer coefficient in a BCR can be defined as,

$$k_L a_i = b_{0,OP1} U_G^{b_{1,OP1}} \quad (5)$$

$U_G$ is the superficial gas velocity, $b_{0,OP1}$ and $b_{1,OP1}$ are parameters that can be affected by the liquid phase properties and type of sparger. Values of 0.467 and 0.82 were respectively proposed by [35] for tap water and salt solutions with $U_G$ ranging between 0.002-0.08 m/s.
**Continuous stirred tank reactor (CSTR).** The volumetric mass transfer coefficient in a CSTR can be defined as,

$$k_L a_i = b_{0,OP2} \left( \frac{Power}{V_{liq}} \right)^{b_{1,OP2}} U_G^{b_{2,OP2}} \quad (6)$$

$b_{0,OP2}$, $b_{1,OP2}$, and $b_{2,OP2}$ are constant parameters. The values of $b_{1,OP2}$ and $b_{2,OP2}$ present variations: $0.4 < b_{1,OP2} < 1$ and $0 < b_{2,OP2} < 0.7$. The value of $b_{0,OP2}$, however, is not reported in most cases, but it is highly correlated to $b_{1,OP2}$ because $(Power/V_{liq})$ range from $1 \times 10^3$ - $1 \times 10^4$ [36]. This means that if small variations are introduced into $b_{1,OP2}$, $b_{0,OP2}$ would show large changes.

The mass transfer coefficient $k_L$ is proportional to the square root of the diffusivity in the liquid phase [37]. In a biomethanation process, some substrates (H₂, CO, and CO₂) are initially present in a gaseous form and must thus transfer from the gas phase to the liquid phase before getting involved in the biochemical reaction. As a result, the biochemical reaction rates might be in the end limited by the mass transfer rate. When comparing different gases, considering both the $k_L a$ and the solubility it appears that the lowest value is that of $k_L a_{H_2}$. For this reason, it can be claimed that if mass transfer becomes the limiting phenomena, it will be the mass transfer rate of H₂ which limit the bioreaction rate [17,38]. Therefore, only $k_L a_{H_2}$ was calculated using Eq. (5) and Eq. (6). The rest were calculated as,

$$k_L a_i = k_L a_{H_2} \sqrt{\frac{D_i}{D_{H_2}}} \quad (7)$$

In this equation the sub-index $i = CH_4, CO, CO_2$. $D_{H_2}$ is the diffusion of H₂ and $D_i$ is the diffusion of CH₄, CO, and CO₂.

An interesting analysis can be performed by increasing $q_{gas}^{in}$ in the BCR, and the Power in the CSTR. For the case of the BCR, it was assumed that the organic loading rate and gas loading rate were the same as proposed in the operational condition one (OP1) (Section 5, Table 3). Fig. 2-A displays the variation of $kLa_i$ with the change of $q_{gas}^{in}$ in a range from $1 \times 10^{-1}$ to $1 \times 10^2 L/d$. All the $kLa_i$ increased with a different rate, i.e., $kLa_{H_2}$ was the higher followed by $kLa_{CO}$, $kLa_{CO_2}$, and $kLa_{CH_4}$. In general, the $U_G$ ranged from $3.2 \times 10^{-3}$ to $3.5 m/d$ for all the $kLa_i$. Values of 23.06 and 39.68 $1/d$ were
obtained for $kLa_{CH4}$ and $kLa_{H2}$ with $U_G$ of 3.4 m/d, while 25.89 1/d was achieved for $kLa_{CO}$ and $kLa_{CO2}$. For the case of the CSTR, the organic loading rate and gas loading rate were assumed to be the same as proposed in operational condition two (OP2) (Section 5, Table 3). The $q_{gas}^{in}$ was fixed at the maximum value proposed for OP2. Fig. 2-B shows the variation of the $kLa_i$ concerning Power which was varied from 1 to $1 \times 10^2$ W. A ($Power/V_{liq}$) between $9.52 \times 10^{-2}$ to 9.61 W/L was obtained from these values. Similar to the BCR, all the $kLa_i$ increased at a different rate, i.e., $kLa_{H2}$ displayed the largest value, followed by $kLa_{CO}$, $kLa_{CO2}$, and $kLa_{CH4}$. Values of 313 and 539 1/d were obtained for $kLa_{CH4}$ and $kLa_{H2}$ in the final stage, while a value of 352 1/d was achieved for $kLa_{CO}$ and $kLa_{CO2}$.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Variation of the volumetric mass transfer coefficient: (A) concerning the added gas flow rate $1 \times 10^{-1} \le q_{gas}^{in} \le 1 \times 10^{-2}$ L/d and $U_G$ for a BCR (Eq. (5)), and (B) respect to the power input $1 \le Power \le 1 \times 10^2$ W for a CSTR (Eq. (6)).</figcaption></figure>

### 3.3. Constitutive Equations

The total ($q_{gas}$) and specific ($q_{gas,i}$) outlet gas flow rates can be calculated with Eqs. (8) and (9) [12].

$$q_{gas} = K_P(P_{gas} - P_{atm}) \qquad (8)$$
$$260 \qquad q_{gas,i} = q_{gas} \left( \frac{P_{gas,i}}{P_{gas}} \right) \qquad (9)$$

In these equations the sub-index $i = CH_4, H_2, CO, CO_2, K_p$ is a parameter related to the friction in the gas outlet, i.e., this parameter fixes the pressure drop at the gas outlet. $P_{atm}$, $P_{gas}$ and $P_{gas,i}$, are the atmospheric, total, and partial pressure of component $i$, respectively. The term $(P_{gas,i}/P_{gas})$ is the molar fraction of component $i$, in which the partial $(P_{gas,i})$, and total pressure $(P_{gas})$ are calculated by using Eqs. (10) – (11).

$$266 \qquad P_{gas,i} = \frac{S_{gas,i}RT}{\gamma_{COD,i}} \qquad (10)$$

$$267 \qquad P_{gas} = \sum P_{gas,i} \qquad (11)$$

where $R$ and $T$ are the Universal gas constant and the temperature of the process, respectively.

For each component $i$ in the gas phase, the percent composition of each $i$ gas ($p_{gas,i}$) was calculated based on the molar fraction, Eq. (12).

$$271 \qquad p_{gas,i} = \left( \frac{P_{gas,i}}{P_{gas}} \right) \cdot 100\% \qquad (12)$$

To model the conversion rates process, the biochemical reaction rates ($\mu_k$) are considered by Monod kinetics for substrate consumption with inhibition, while biomass decay ($\mu_{k,dec}$) follows first-order kinetics. This allows the representation of the cell growth and death associated with each biomass involved in the process. Their mathematical representation is expressed as,

$$276 \qquad \mu_k = \frac{\mu_{m,k} S_{liq,j}}{KS_k + S_{liq,j}} X_k I_{H2,k} I_{CO,k} \qquad (13)$$

$$277 \qquad \mu_{k,dec} = K_{k,dec} X_k \qquad (14)$$

where $\mu_{m,k}$ is the maximum specific growth rate, $KS_k$ the saturation constant, and $K_{k,dec}$ the decay biomass constant. The term $I_{H2,k}$ in Eq. (13) concerns inhibition with respect to $H_2$, which affects the reaction rates of butyrate, propionate, and acetate. The term $I_{CO,k}$ concerns inhibition with respect to CO, which occurs for the reaction rates regarding acetate and $H_2$. Eqs (15) – (16) represent these inhibitions.
$$283 \qquad I_{H2,k} = \frac{1}{1 + S_{liq,H2}/KI_{H2,k}} \qquad (15)$$

$$284 \qquad I_{CO,k} = \frac{1}{1 + S_{liq,CO}/KI_{CO,k}} \qquad (16)$$

$KI_{H2,k}$ are the inhibition constants over butyrate, propionate, and acetate due to $H_2$, meanwhile $KI_{CO,k}$ are the inhibition constants of the effect of CO over acetate and $H_2$.

## 4. Global Sensitivity Analysis Techniques

Biological kinetic models are constituted of several parameters, e.g., stoichiometric, biochemical, and physicochemical parameters, which present a certain degree of uncertainty. Sensitivity analysis (SA) is a powerful alternative to determine how the uncertainty of the model inputs or parameters influence the outputs [39–42]. Sensitivity Analysis was developed to identify the contribution of each uncertainty of the inputs on the outputs [43–45]. Techniques for SA can be classified into local and global [46,47]. In the context of this work, we will focus on Global sensitivity analysis (GSA).

In GSA, the model is studied from a probabilistic point of view. The exploration of the entire range of variation of the model parameters is considered using a probability density function associated with each input parameter and repeated simulations of the model [42,47,48]. GSA could employ regression methods, screening methods, and variance-based methods [49].

### 4.1. Sobol's Method

The Sobol method [50] is an interesting variance-based method in which the variance of the model output can be decomposed into partial variances that represent the contribution of the inputs over the overall uncertainty of the model output [40,42,46,47,50].

Consider the model define by $\xi$, Eq. (17).

$$303 \qquad Y = [S_{liq,j}, X_k, S_{gas,i}] = \xi(\theta) \qquad (17)$$

where $Y \in \mathbb{R}^m$ is the model output of interest, and $\theta \in \mathbb{R}^n$ is a $n$-dimensional parameter vector defined as $\theta = (\theta_1, \theta_2, ..., \theta_n)$ and characterized by a probability density function (PDF).

The function $\xi(\theta)$ can be decomposed into summands of different dimensions, Eq. (18).
$$307 \qquad Y = \xi_0 + \sum_{i=1}^{n} \xi_i(\theta_i) + \sum_{1 \le i < j \le n} \xi_{ij}(\theta_i, \theta_j) + \dots + \xi_{1\dots n}(\theta_i, \dots, \theta_n) = \sum_{\mathbf{u} \subseteq \{1, \dots, p\}} \xi_{\mathbf{u}}(\theta_{\mathbf{u}}) \qquad (18)$$

where:

$$308 \qquad \left\{ \begin{array}{l} \xi_0 = \mathbb{E}[Y] \\ \xi_i(\theta_i) = \mathbb{E}[Y|\theta_i] - \xi_0 \\ \xi_{ij}(\theta_i, \theta_j) = \mathbb{E}[Y|\theta_i, \theta_j] - \xi_0 - \xi_i - \xi_j \end{array} \right. \qquad (19)$$

$\xi_0$ is the mean of the function, $\xi_i(\theta_i)$ and $\xi_{ij}(\theta_i, \theta_j)$ are the expectation terms of increasing order and the conditional expectations defined recursively. This decomposition is unique provided that the inputs are independent and the individual terms are square integrable. Squaring Eq. (18) and integrating, we can get the so-called ANOVA decomposition as,

$$310 \qquad V(Y) = \sum_{\mathbf{u}} V(\xi_{\mathbf{u}}(\theta_{\mathbf{u}})) \quad \text{for } \mathbf{u} \subset \{1, \dots, n\} \qquad (20)$$

where $V(\xi_{\mathbf{u}}(\theta_{\mathbf{u}}))$ express the conditional variance for the subvector $\theta_{\mathbf{u}}$, containing the variables whose indices are indicated by the subset $\mathbf{u}$.

The variance of the output can be decomposed into terms depending on the parameters and their interactions. The variance-based sensitivity index (SI) associated to the subset $\mathbf{u}$ is defined as the ratio between the contribution given by the interaction among the components of $\mathbf{u}$ for the model variance, and the total variance, Eq. (21).

$$321 \qquad SI_{\mathbf{u}} = \frac{V(\xi_{\mathbf{u}}(\theta_{\mathbf{u}}))}{V(Y)} \qquad (21)$$

Based on this, for $\mathbf{u} \subset \{1, \dots, n\}$, and $\mathbf{u} \neq \mathbf{0}$,

$$322 \qquad \sum_{\mathbf{u}} SI_{\mathbf{u}} = \sum_{i=1}^{n} SI_i + \sum_{1 \le i < j \le n} SI_{ij} + \dots + SI_{1\dots n} = 1 \qquad (22)$$

The term $SI_i$ is the first-order sensitivity index, which measures the fraction of the total output variance explained by the parameter $\theta_i$ alone as,

$$326 \qquad SI_i = \frac{V(\xi_i(\theta_i))}{V(Y)} \quad i = 1, \dots, n \qquad (23)$$
Similarly, $SI_{ij}$ is the second order-sensitivity index that measures the amount of variance caused by the interaction between the parameters $\theta_i$ and $\theta_j$ as,

$$ SI_{ij} = \frac{V(\xi_{ij}(\theta_{ij}))}{V(Y)} \quad 1 \le i < j \le n \qquad (24) $$

It is possible to construct the $SI$ for all orders until the $n^{th}$ order index $SI_{1,...,n}$, which represents the contribution of the interactions between all the parameters in $\theta$. To measure the full contribution of the $i^{th}$ random parameter $\theta_i$ for the total variance either by its single effect or by its interaction with others, we use the total Sobol' indices,

$$ SI_i^T = \sum_{\substack{u \in \{1,...,n\} \\ i \in u}} SI_{ij} + \dots + SI_{1...n} \quad i = 1, ..., n \qquad (25) $$

Eq. (25) indicates that the total sensitivity index does not only include the marginal contribution of $\theta_i$ to the variance of the output, but it also contains its cooperative contribution with all the other inputs. The total sensitivity index removes the parameter $\theta_i$ from the analysis and allocates the resulting variance reduction to that parameter.

## 4.2. Morris Method

The Morris method is the most well-known [51] screening method which performs SA by analyzing one-factor-at-a-time (OAT). This is generally used when the number of model parameters is large and the computation of model simulations is expensive. This type of method provides qualitative sensitivity measures, ranking the factors in accordance with their importance. Nevertheless, it does not quantify the importance of one factor with respect to another [52]. The Morris method applied to parameter sensitivity discretizes the space of each parameter and performs a given number of OAT designs. These designs and variation directions are randomly chosen from the parameter space. The repetition of these steps allows the estimation of elementary effects (EE) for each input, which represents the relative difference between the outputs and the $i^{th}$ parameter disturbance [43,46,48,51,53].

Consider a trajectory in the parameter space as,

$$ \theta_i = \theta_{i-1} + e_i \Delta_i \qquad (26) $$
According to Eq. (26), the elementary effects can be calculated as follows,

$$EE_i = \frac{\xi(\theta_{i-1} + e_i\Delta_i) - \xi(\theta_i)}{\Delta_i} \qquad (27)$$

where $\theta_i$ is a sample of input $\theta$ and $\xi(\theta_i)$ is the corresponding model output. $e_i$ is a vector of zeros but with a one in its $i^{th}$ component, $\Delta_i$ is a step between two consecutive input space points of the trajectory. The term $(\theta_{i-1} + e_i\Delta_i)$ represents a new sample by moving the $i^{th}$ input from $\theta_i$ to $\theta_i + \Delta_i$, with the respective model output $\xi(\theta_i + e_i\Delta_i)$.

The index of $EE_i$ expresses the ratio of the change of Y (output) when the $i^{th}$ parameter $\theta_i$ is given a certain change $\Delta_i$. Then, $EE_i$ can measure the effect of $\theta_i$ in a given scope of output Y. Once N elementary effects per input are available ($EE_i^j$ ($i = 1, 2, ..., n, j = 1, 2, ..., N$)), some statistics of the distribution of the elementary effects can be chosen to indicate the global sensitivity of parameter $\theta_i$. The sensitivity measures expressed in terms of means $\mu_i^*$, and standard deviations $\sigma_i$ are defined as Eqs. (28) – (29).

$$\mu_i^* = \frac{\sum_{j=1}^{r} |EE_i^j|}{N} \qquad (28)$$

$$\sigma_i = \sqrt{\frac{\sum_{j=1}^{r} (EE_i^j - \mu_i)^2}{N - 1}} \qquad (29)$$

The sensitivity measures $\mu_i$ and $\sigma_i$ are the mean and the standard deviation of the distribution of the elementary effects, respectively. $\mu_i^*$ measures the influence of the $j^{th}$ parameter on the output. $\sigma_i$ is a measure of non-linear and/or interaction effects of the $j^{th}$ parameter. $EE_i$ measures the effect of $\theta_i$ on the output at a local point, the mean $\mu_i^*$ can be used to assess the overall effect of $X_i$ on the output. A high value of $\mu_i^*$ indicates that the parameter $\theta_i$ has a more important effect on the output. A high value of $\sigma_i$ indicates that the elementary effect of $\theta_i$ varies significantly from one to another, which shows that the value of $EE$ is strongly influenced by the selected sample points.

## 5. Bioreactors Operating Conditions

Two experimental data sets from the literature were used in this study to derive the model extension.

The first data set was taken from [16] where the experiments were developed at OP1. The second data
set was obtained from [54] whose experiment was carried out at OP2. Both operating conditions are reported in Table 2.

The measured outputs were different for each operating condition. For OP1, these corresponded to the outlet flow rates of CH₄, H₂, and CO; while for OP2 the measured outputs were the percent of CH₄, H₂, CO, and CO₂ in the gas phase.

<table><caption>**Table 2.** Operational conditions from OP1 [16] and OP2 [54].</caption><thead><tr><th>Operational conditions</th><th>OP1</th><th>OP2</th></tr></thead><tbody><tr><td>Reactor type</td><td>BCR</td><td>CSTR</td></tr><tr><td>Temperature (°C)</td><td>37</td><td>55</td></tr><tr><td>Working volume (L)</td><td>37.5</td><td>10.5</td></tr><tr><td>HRT (d)</td><td>20</td><td>21</td></tr><tr><td>Experimental time (d)</td><td>207</td><td>150</td></tr><tr><td>Inlet liquid flow rate (L/d)</td><td>1.9</td><td>0.5</td></tr></tbody></table>

HRT: hydraulic retention time.

In OP1, the loading consisted of two additions: glucose and syngas. The organic loading rate (OLR) of glucose was kept at 0.5 $gGlu/L_r/d$, where $L_r$ represents the volume of the reactor. Syngas containing 50% v/v of H₂ and CO ($H_2/CO \approx 1$) was added into the reactor with a continuous flow, but at different rates after a first reference stage. For OP2, the liquid fraction of the OLR was a mixture of primary sludge and activated ticked-disintegrated waste (volume ratio of 3:1) which was fed with a flow rate varying in time. Syngas was also added into the reactor at diverse rates after a reference stage. In this case, however, the syngas contained 55% v/v of H₂ and 45% v/v CO ($H_2/CO \approx 0.55/0.45$). Table 3 reports the stages in which the substrates were added to the reactors.

To standardize the units for both operating conditions, the OLR is expressed as $gCOD/L_r/d$ and the gas loading rate is expressed as $L/L_r/d$. In OP1, the OLR was expressed in $gCOD/L_r/d$ by multiplying the value 0.5 $gGlu/L_r/d$, by a factor of 1.07 $gCOD/gGlu$. In OP2, a mixture of different substrates was used. Those mixtures are normally expressed in terms of volatile solids (VS). Therefore, an equivalence between COD and VS must be estimated. Some authors mentioned 1.42 $gCOD/gVS$
for activated sludge [55], while others indicated values between 1.6 – 1.7 gCOD/gVS [56]. In this case, an estimation using the function *fmincon* from MATLAB® was performed by minimizing the RMSE with a lower and upper value of 1.42 and 1.7 gCOD/gVS, respectively. A value of 1.62 gCOD/gVS was obtained. Additionally, a value of 0.025 mol/gCOD was used for the CO₂ balance (lower than the one proposed for sugar 0.0313 mol/gCOD). To model the mass transfer for OP2, a value of 0.69 was used for *b*₁,OP1 in Eq. (5) since the BCR structure is not known in detail.

<table><caption>**Table 3.** Syngas flow rate, gas loading rate, and organic loading rate from OP1 [16] and OP2 [54].</caption><thead><tr><th>Stage</th><th>Time (d)</th><th>Gas flow rate (L/d)</th><th>Gas loading rate (L/L<sub>r</sub>/d)</th><th>Organic loading rate (gCOD/L<sub>r</sub>/d)</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">OP1</td></tr><tr><td>Reference</td><td>1-32</td><td>0.0</td><td>0.0</td><td></td></tr><tr><td>I</td><td>33-64</td><td>7.5</td><td>0.2</td><td></td></tr><tr><td>II</td><td>65-101</td><td>7.5</td><td>0.2</td><td rowspan="3">0.53</td></tr><tr><td>III</td><td>102-135</td><td>15.0</td><td>0.4</td></tr><tr><td>IV</td><td>136-171</td><td>37.5</td><td>1.0</td></tr><tr><td>V</td><td>172-207</td><td>37.5</td><td>1.0</td><td></td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">OP2</td></tr><tr><td>Reference</td><td>1-36</td><td>0.0</td><td>0.0</td><td>3.08</td></tr><tr><td>I</td><td>36-51</td><td>3.15</td><td>0.3</td><td>3.72</td></tr><tr><td>II</td><td>51-81</td><td>7.35</td><td>0.7</td><td>3.24</td></tr><tr><td>III</td><td>81-118</td><td>10.5</td><td>1.0</td><td>3.09</td></tr><tr><td>IV</td><td>118-130</td><td>15.75</td><td>1.5</td><td>2.84</td></tr></tbody></table>

## 6. Model Calibration and Validation

In this study, model simulations were implemented in MATLAB® and run using a computer with Intel® Core i7 8665U 2.11 GHz and 16 GB RAM. The ADM1_ME was calibrated using the measured
outputs mentioned in the previous section, CH₄, H₂, and CO outlet flow rates for OP1, and CH₄ H₂, CO, and CO₂ percent in the gas phase for OP2.

Sobol's method was implemented using the toolbox: Global sensitivity and uncertainty analysis (GSUA) [57], whereas the Morris method was performed by the toolbox: Sensitivity analysis-Morris method (advanced) [58]. The two abovementioned SA methods were performed to identify the effects of a change in the parameters on the model outputs.

The ADM1_ME is described by 60 parameters comprising stoichiometric and kinetic parameters (see supplementary material Table S.1). Among those parameters, a group of 26 parameters was selected to analyze their impact on the gas flow rates for OP1 and gases percent for OP2 (see supplementary material Table S.1). Those comprised mostly kinetic parameters and the parameters where the estimation is more uncertain. The Morris method performs analysis upon several simulation runs, e.g., $(n + 1) \cdot 10$. In this case, however, a larger number of simulation runs were executed $(n + 1) \cdot 100$, i.e., 2700 to guarantee a good sampling in the distribution of the parameter domain. Since the Sobol' method is computationally expensive with respect to the Morris method, only 200 simulation runs were evaluated per output.

The toolbox GSUA for Sobol' method allowed the computation of the first-order sensitivity indices with a scalar characteristic ($SI_{is}$) as,

$$SI_{is} = \xi(SSE) = \xi(\sum(Y_m(\theta_i) - Y_m(\theta_{i,nom}))^2) \qquad (30)$$

where the sum squared error (SSE) was the scalar characteristic, measured between the output variable ($Y_m(\theta_i)$) calculated with the varied parameter $\theta_i$, sampled with a uniform distribution, and the output variable ($Y_m(\theta_{i,nom})$) calculated with the nominal value of the parameter $\theta_{i,nom}$.

The most sensitive parameters founded with SA were selected to recalibrate the model, while the rest of the parameters were fixed to the nominal values reported in the ADM1 [11].

Parameter estimation was performed to minimize the adapted root mean square error (RMSE$_{est}$) reported in Eq. (31). The outputs for OP1 were the outlet gas flow rates, $q_{gas,CH4}$, $q_{gas,H2}$, and $q_{gas,CO}$. Nevertheless, for OP2 the outputs were the gases percent, $p_{gas,CH4}$, $p_{gas,H2}$, $p_{gas,CO}$, and $p_{gas,CO2}$. Therefore, Eq. (31) considers $q_{gas,i}$ for OP1 and $p_{gas,i}$ for OP2.
$$434 \quad RMSE_{est} = \sqrt{\frac{\omega_{OP1}}{n_{OP1}} \sum_{i=1}^{3} \left( \frac{|\hat{q}_{gas,i} - q_{gas,i}|}{\max(q_{gas,i})} \right)^2 + \frac{\omega_{OP2}}{n_{OP2}} \sum_{i=1}^{4} \left( \frac{|\hat{p}_{gas,i} - p_{gas,i}|}{\max(p_{gas,i})} \right)^2} \quad (31)$$

Sub-index *i* denotes CH₄, H₂, CO, and CO₂. Note that only the first three are used for OP1. *n*<sub>OP1</sub> and *n*<sub>OP2</sub> are the number of observations in OP1 and OP2. *ω*<sub>OP1</sub> and *ω*<sub>OP2</sub> are the weights to trade-off the estimation of OP1 and OP2.

Model validation was carried out with the results from OP1 (stages IV and V) and the results from OP2 (stages III and IV). The coefficient of determination (*R*²), and the root mean squared error (*RMSE*) were used as criteria to qualify parameter estimation, Eqs. (32) – (33).

$$441 \quad R_i^2 = 1 - \frac{\sum (y_{exp} - y_{mod})^2}{\sum (y_{exp} - \bar{y}_{exp})^2} \quad (32)$$

$$442 \quad RMSE_i = \sqrt{\frac{1}{n_{data set}} \sum (y_{exp} - y_{mod})^2} \quad (33)$$

In the former equations, *i* denotes CH₄, H₂, and CO for OP1, and denotes CH₄, H₂, CO, and CO₂ for OP2. *n*<sub>data set</sub> is the number of observations for the operating condition, *y*<sub>exp</sub>, *ȳ*<sub>exp</sub>, and *y*<sub>mod</sub> are the experimental, mean, and model data, respectively.

Confidence intervals were determined for each estimated parameter by computing the Global Sensitivity Information Matrix (*GSIM*). This matrix is based on the first-order sensitivity indices (*SI*<sub>*i*</sub>) of the Sobol' method and is analogous to the Fisher Information Matrix (*FIM*) [59,60]. The *GSIM* is calculated as,

$$450 \quad GSIM = \sum_{t=1}^{n_{data set}} [Q^T(t_t) W_t^{-1} Q(t_t)] \quad (34)$$

where $W_i^{-1}$ is a weighting matrix usually chosen as the measurement error covariance matrix, and $Q(t_t)$ is defined as,

$$453 \quad Q(t_t) = \begin{bmatrix} SI_{1,1} & SI_{1,2} & \dots & SI_{1,p} \\ SI_{2,1} & SI_{2,2} & \dots & SI_{2,p} \\ \vdots & \vdots & \ddots & \vdots \\ SI_{i,1} & SI_{i,2} & \dots & SI_{i,p} \end{bmatrix} \quad (35)$$
In this case, $SI_{i,p}$ was used to indicate the $SI_i$ for each model output $p$. The variance of each parameter $\theta_i$ can be approximated by $\sigma^2(\theta_i) \approx GSIM_{ii}^{-1}$ and used to evaluate the 95% confidence intervals as: $\theta_i \pm 1.96 \cdot \sigma(\theta_i)$ where $\sigma(\theta_i)$ is the standard deviation.

# 7. Results and Discussion

## 7.1. Sensitivity Analysis

The first-order sensitivity indices of Sobol' Method were computed for the 26 selected parameters and evaluated for each output and operating condition. Table 4 summarizes the results. The sum of the variance is considered to be around 100%. This normalization helps to identify which parameters contribute the most to the total variance, i.e., the most sensitive parameters (see supplementary material Fig. S.1 - S.2).

<table><caption>**Table 4.** First-order sensitivity index with a scalar characteristic ($SI_{is}$) with the Sobol' method from OP1 [16] and OP2 [54].</caption><thead><tr><th rowspan="2">Parameter / Model outputs</th><th colspan="4">SI<sub>is</sub> for OP1</th><th colspan="4">SI<sub>is</sub> for OP2</th></tr><tr><th>q<sub>gas,CH4</sub></th><th>q<sub>gas,H2</sub></th><th>q<sub>gas,CO</sub></th><th>q<sub>gas,CO2</sub></th><th>p<sub>gas,CH4</sub></th><th>p<sub>gas,H2</sub></th><th>p<sub>gas,CO</sub></th><th>p<sub>gas,CO2</sub></th></tr></thead><tbody><tr><td>Y<sub>su</sub></td><td>2.8</td><td>3.0</td><td>6.2</td><td>2.5</td><td>0.7</td><td>8.4</td><td>3.6</td><td>4.1</td></tr><tr><td>Y<sub>bu</sub></td><td>4.1</td><td>5.9</td><td>6.7</td><td>8.0</td><td>1.3</td><td>2.9</td><td>2.5</td><td>1.7</td></tr><tr><td>Y<sub>pro</sub></td><td>5.6</td><td>6.4</td><td>3.5</td><td>2.6</td><td>1.4</td><td>1.5</td><td>3.5</td><td>3.2</td></tr><tr><td>Y<sub>ac</sub></td><td>4.4</td><td>7.2</td><td>1.8</td><td>2.4</td><td>1.4</td><td>5.3</td><td>4.7</td><td>2.2</td></tr><tr><td>Y<sub>CO</sub></td><td>3.0</td><td>2.5</td><td>5.7</td><td>2.9</td><td>5.1</td><td>1.7</td><td>5.0</td><td>1.4</td></tr><tr><td>Y<sub>H2</sub></td><td>3.9</td><td>4.7</td><td>1.8</td><td>3.1</td><td>8.7</td><td>2.7</td><td>2.8</td><td>2.5</td></tr><tr><td>μ<sub>m,su</sub></td><td>3.1</td><td>9.0</td><td>2.4</td><td>10.0</td><td>3.3</td><td>2.7</td><td>5.7</td><td>8.4</td></tr><tr><td>μ<sub>m,bu</sub></td><td>3.2</td><td>3.6</td><td>3.9</td><td>9.6</td><td>0.7</td><td>3.1</td><td>2.1</td><td>2.5</td></tr><tr><td>μ<sub>m,pro</sub></td><td>3.0</td><td>3.4</td><td>3.1</td><td>1.7</td><td>5.3</td><td>11.8</td><td>8.8</td><td>0.3</td></tr><tr><td>μ<sub>m,ac</sub></td><td>4.0</td><td>2.2</td><td>3.0</td><td>4.4</td><td>1.2</td><td>2.5</td><td>2.8</td><td>6.0</td></tr><tr><td>μ<sub>m,CO</sub></td><td>3.1</td><td>4.0</td><td>3.0</td><td>2.7</td><td>1.2</td><td>4.4</td><td>11.5</td><td>3.1</td></tr></tbody></table>
<table><tr><td>μ<sub>m,H2</sub></td><td>3.3</td><td>5.7</td><td>4.7</td><td>2.7</td><td>0.6</td><td>3.2</td><td>4.0</td><td>11.3</td></tr><tr><td>f<sub>ac,CO</sub></td><td>3.7</td><td>2.5</td><td>4.3</td><td>3.7</td><td>1.3</td><td>1.8</td><td>1.3</td><td>5.4</td></tr><tr><td>KI<sub>H2, ac</sub></td><td>4.0</td><td>3.3</td><td>2.6</td><td>2.1</td><td>44.2</td><td>2.0</td><td>3.9</td><td>14.6</td></tr><tr><td>KI<sub>H2, bu</sub></td><td>4.9</td><td>3.9</td><td>3.3</td><td>6.7</td><td>0.8</td><td>3.4</td><td>4.8</td><td>0.4</td></tr><tr><td>KI<sub>H2, pro</sub></td><td>3.9</td><td>3.3</td><td>3.7</td><td>6.2</td><td>1.3</td><td>7.0</td><td>3.3</td><td>2.7</td></tr><tr><td>KI<sub>CO, ac</sub></td><td>4.4</td><td>3.0</td><td>4.5</td><td>2.5</td><td>3.3</td><td>2.5</td><td>3.8</td><td>3.6</td></tr><tr><td>KI<sub>CO, H2</sub></td><td>3.8</td><td>1.5</td><td>5.2</td><td>3.1</td><td>0.4</td><td>1.4</td><td>3.3</td><td>3.0</td></tr><tr><td>KS<sub>su</sub></td><td>3.5</td><td>2.6</td><td>4.0</td><td>3.4</td><td>1.7</td><td>8.6</td><td>2.5</td><td>2.7</td></tr><tr><td>KS<sub>bu</sub></td><td>4.2</td><td>1.7</td><td>3.3</td><td>2.1</td><td>2.1</td><td>3.6</td><td>2.0</td><td>0.9</td></tr><tr><td>KS<sub>pro</sub></td><td>4.7</td><td>4.9</td><td>3.4</td><td>2.3</td><td>0.8</td><td>2.2</td><td>3.1</td><td>2.2</td></tr><tr><td>KS<sub>ac</sub></td><td>4.9</td><td>2.4</td><td>3.0</td><td>2.8</td><td>0.7</td><td>2.1</td><td>1.6</td><td>2.7</td></tr><tr><td>KS<sub>CO</sub></td><td>2.4</td><td>2.4</td><td>3.9</td><td>3.9</td><td>4.8</td><td>3.1</td><td>3.2</td><td>2.7</td></tr><tr><td>KS<sub>H2</sub></td><td>3.0</td><td>2.9</td><td>5.1</td><td>3.0</td><td>1.8</td><td>3.2</td><td>3.7</td><td>5.4</td></tr><tr><td>k<sub>L</sub>a<sub>H2</sub></td><td>4.8</td><td>5.9</td><td>2.6</td><td>2.9</td><td>3.6</td><td>1.5</td><td>3.7</td><td>5.4</td></tr><tr><td>K<sub>P</sub></td><td>4.2</td><td>2.3</td><td>5.1</td><td>2.8</td><td>2.3</td><td>7.5</td><td>2.9</td><td>1.5</td></tr><tr><td>Σ</td><td>99.9</td><td>100.2</td><td>99.8</td><td>100.1</td><td>100.0</td><td>100.1</td><td>100.1</td><td>99.9</td></tr></table>

466

Sobol' method allowed the determination of the first-order sensitivity index given the scalar characteristic. A threshold value of 5% was proposed to consider which parameters were sensitive.

467

468

Concerning OP1, a value of 5.6 was obtained for the first-order sensitivity index of $Y_{pro}$ with respect to $q_{gas,CH4}$. For $q_{gas,H2}$, the parameters $μ_{m,su}$, $Y_{ac}$, $Y_{pro}$, $Y_{bu}$, $k_L a_{H2}$, and $μ_{m,H2}$ were found to be the most sensitive. Regarding $q_{gas,CO}$, the parameters reporting the highest values were $Y_{bu}$, $Y_{su}$, $Y_{CO}$, $KI_{CO,H2}$, $KS_{H2}$, and $K_P$. Finally, the parameters $μ_{m,su}$, $μ_{m,bu}$, $Y_{bu}$, $KI_{H2, bu}$, and $KI_{H2, pro}$ were found to be the most sensitive for $q_{gas,CO2}$.

469

470

471

472

473

For OP2, the parameters $KI_{H2,ac}$, $μ_{m,pro}$, $Y_{CO}$, and $Y_{H2}$ were reported as the most sensitivities for the $p_{gas,CH4}$. Regarding $p_{gas,H2}$, the parameters $μ_{m,pro}$, $KS_{su}$, $Y_{su}$, $K_P$, $KI_{H2, pro}$, and $Y_{ac}$ were found to be

474

475
the most sensitive, while the parameters reporting the highest values for $p_{gas,CO}$ were $\mu_{m,CO}$,
$\mu_{m,pro}$, $K_{S_{su}}$, $\mu_{m,su}$, and $Y_{CO}$. Concerning $q_{gas,CO_2}$, the most sensitive parameters were $KI_{H_2,ac}$,
$\mu_{m,H_2}$, $\mu_{m,su}$, $\mu_{m,ac}$, $\beta_{ac,CO}$, $K_{S_{H_2}}$, and $k_L a_{H_2}$.

In OP2, $K_{S_{su}}$ presented an effect on the model outputs $q_{gas,CH_4}$ and $q_{gas,CO_2}$, while $KI_{H_2,ac}$ affected $q_{gas,H_2}$ and $q_{gas,CO}$. The growth rate $\mu_{m,pro}$ presented an effect over $q_{gas,CO}$ and $q_{gas,CO_2}$, while $\mu_{m,su}$ did it over $q_{gas,CH_4}$ and $q_{gas,H_2}$. The model outputs $q_{gas,H_2}$ and $q_{gas,CO_2}$ were also affected by $Y_{CO}$.

The other high-influence parameters reported an effect in just one output. Therefore, the selection of unique candidates for parameter estimation was not straightforward. Nevertheless, parameters that appeared as sensitive for most of the model outputs were considered good candidates for parameter estimation, e.g., $Y_{bu}$, $Y_{pro}$, $Y_{CO}$, $\mu_{m,su}$, $\mu_{m,pro}$, $KI_{H_2,ac}$, $K_{S_{su}}$, and $K_p$.

The elementary effects of the Morris method were computed for all 26 parameters. Figs. 3 and 4 show only the high influence factors identified for each output ordered by ascending maximum (see supplementary material Table S.2 for more details). For OP1, 12, 10, 8, and 5 parameters were determined as high-influence factors over $q_{gas,CH_4}$, $q_{gas,H_2}$, $q_{gas,CO}$, and $q_{gas,CO_2}$, respectively. Differently, for OP2, 10, 11, 7, and 18 parameters were determined as high influence factors over $p_{gas,CH_4}$, $p_{gas,H_2}$, $p_{gas,CO}$, and $p_{gas,CO_2}$, respectively.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Morris sensitivity analysis with OP1 over the model outputs: (A) $q_{gas,CH_4}$, (B) $q_{gas,H_2}$, (C) $q_{gas,CO}$, and (D) $q_{gas,CO_2}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Morris sensitivity analysis with OP2 over the model outputs: (A) $p_{gas,CH4}$, (B) $p_{gas,H2}$, (C) $p_{gas,CO}$, and (D) $p_{gas,CO2}$.</figcaption></figure>

As both operating conditions represented different output variables (OP1: outlet gas flow rates and OP2: gases percent), the sensitivity with respect to each factor could be different. Additionally, the operating conditions are not the same. OP1 was carried out under mesophilic temperatures while OP2 was developed under thermophilic conditions. These conditions could affect the kinetics at which microorganisms are converting the substrates. Therefore, the best trade-off was pursued with both sets of parameters. In OP1 parameters, such as $\mu_{m,su}$, $Y_{su}$, $Y_{ac}$, $\mu_{m,ac}$, and $K_{s_{su}}$ presented an influence over all the model outputs, and $Y_{pro}$, $K_p$, and $\mu_{m,pro}$ affected the model outputs $q_{gas,CH4}$, $q_{gas,H2}$, and $q_{gas,CO}$. $K_{s_{ac}}$ and $K_{s_{pro}}$ only influenced $q_{gas,CH4}$ and $q_{gas,H2}$. The other high influence parameters presented an effect just in one output. In OP2, $\mu_{m,su}$, $K_{s_{H2}}$, $\mu_{m,H2}$, $Y_{su}$ influence all the model outputs. $Y_{ac}$, $\mu_{m,ac}$, and $Y_{H2}$ reported an influence on $q_{gas,CH4}$, $q_{gas,H2}$, and $q_{gas,CO2}$, and $K_{s_{CO}}$, $\mu_{m,CO}$, and
$Y_{CO}$ affected $q_{gas,H2}$, $q_{gas,CO}$, and $q_{gas,CO2}$. $KI_{H2,ac}$, and $KS_{ac}$ presented an influence on $q_{gas,CH4}$ and $q_{gas,CO2}$. The other high influence parameters presented an effect just in one output. Therefore, good candidates for parameter estimation were $Y_{su}$, $Y_{pro}$, $Y_{ac}$, $\mu_{m,su}$, $\mu_{m,pro}$, $\mu_{m,ac}$, $\mu_{m,H2}$, $KS_{su}$, $KS_{ac}$, $KS_{H2}$, and $K_P$.

Both methods were similar in terms of how to measure the parameters contribution to each model output, either by the first-order sensitivity index or by the elementary effects. However, Sobol' method was more computationally expensive. For instance, more than 1h to run 200 simulations per output, while the Morris method ran 2700 simulations per output in less than 1h. Comparing the two sensitivity methods, it was possible to select a group of parameters to be estimated. Table 5 reports the 14 most sensitive parameters of ADM1_ME that affect the model outputs of both operating conditions.

## 7.2. Parameter Estimation and Model Validation

The estimation of the 14 most sensitive parameters was performed using the *fmincon* function from MATLAB® to solve an optimization problem whose objective function was given by Eq. (31). The weights $\omega_{DS1}$ and $\omega_{DS2}$ were manually adjusted to values of 1 and $1 \times 10^3$. To guarantee that the optimum was a local optimum, simulations were run 10 times adding 10% noise to each obtained parameter after the third iteration. The *RMSE*$_{est}$ was minimized to a value of 10.67.

Table 5 shows the estimated parameters and their respective confidence intervals. Parameters such as $Y_{su}$, $Y_{bu}$, $Y_{pro}$, $Y_{ac}$, $Y_{CO}$, $\mu_{m,H2}$, $KS_{su}$, $KS_{ac}$ and $K_P$ have confidence intervals according to their magnitude, and contain the initial value used in the estimation. Parameters such as $\mu_{m,su}$, $KI_{H2,ac}$, and $KS_{H2}$ have smaller confidence intervals, indicating that their values cannot change significantly, i.e., there is a 95% probability of finding their precise value in the given interval. The confidence intervals for $\mu_{m,ac}$ and $\mu_{m,pro}$ are comparable with the mentioned work [15] and the values reported in the original ADM1 [11].

<table><caption>**Table 5.** Estimated parameters for the ADM1_ME.</caption><thead><tr><td>Parameter</td><td>Estimated value</td><td>Confidence Interval</td></tr></thead></table>
<table><tr><td>Y<sub>su</sub></td><td>8.14×10<sup>-2</sup></td><td>±2.9×10<sup>-2</sup></td></tr><tr><td>Y<sub>bu</sub></td><td>6.05×10<sup>-2</sup></td><td>±2.6×10<sup>-2</sup></td></tr><tr><td>Y<sub>pro</sub></td><td>2.81×10<sup>-2</sup></td><td>±1.2×10<sup>-2</sup></td></tr><tr><td>Y<sub>ac</sub></td><td>4.29×10<sup>-2</sup></td><td>±1.8×10<sup>-2</sup></td></tr><tr><td>Y<sub>CO</sub></td><td>2.26×10<sup>-2</sup></td><td>±5.6×10<sup>-3</sup></td></tr><tr><td>μ<sub>m,su</sub></td><td>31.59</td><td>±2.0×10<sup>-1</sup></td></tr><tr><td>μ<sub>m,pro</sub></td><td>10.40</td><td>±1.74</td></tr><tr><td>μ<sub>m,ac</sub></td><td>8.79</td><td>±2.38</td></tr><tr><td>μ<sub>m,H2</sub></td><td>109.47</td><td>±15.09</td></tr><tr><td>KI<sub>H2, ac</sub></td><td>9.75×10<sup>-7</sup></td><td>±2.15×10<sup>-8</sup></td></tr><tr><td>K<sub>s,su</sub></td><td>2.11×10<sup>-2</sup></td><td>±1.5×10<sup>-2</sup></td></tr><tr><td>K<sub>s,ac</sub></td><td>4.96×10<sup>-2</sup></td><td>±1.1×10<sup>-2</sup></td></tr><tr><td>K<sub>s,H2</sub></td><td>1.0270×10<sup>-6</sup></td><td>±7.15×10<sup>-7</sup></td></tr><tr><td>K<sub>p</sub></td><td>4.99×10<sup>4</sup></td><td>±1.8×10<sup>4</sup></td></tr></table>

Table 6 reports the statistical evaluation of the calibration of the ADM1_ME with the two operating conditions. Values higher than 0.90 and 0.91 were determined for $R^2$ with OP1 and OP2. The variables $q_{gas,CH4}$ and $p_{gas,CH4}$ reported the highest values corresponding to 0.95 and 0.98, respectively. In general, the ADM1_ME showed a good fitting to the results of both operating conditions. Concerning the RMSE, values lower than 0.38 and 2.52 were obtained with OP1 and OP2 where the variables $q_{gas,CH4}$, $q_{gas,CO}$, $p_{gas,CH4}$ and $p_{gas,CH4}$ displayed the best results. The model showed better accuracy on these variables, i.e., the variation in error when the model results and the operating conditions are compared is lower for these variables.

The difference in the magnitude between the RMSE for OP1 and OP2 is related to the magnitude of $\Delta q_{gas,i}$ and $\Delta p_{gas,i}$ (numerators in Eq. (31)), and the use of the weights $\omega_{OP1}$ and $\omega_{OP2}$ previously proposed ($\omega_{OP2} = 1 \times 10^3 \omega_{OP1}$) helps to counterbalance the differences in variable magnitude when the optimizer is applied.
<table><caption>**Table 6. Statistical analysis for ADM1_ME calibration with OP1 and OP2.**</caption><thead><tr><th rowspan="2">Criteria</th><th colspan="4">OP1</th><th colspan="4">OP2</th></tr><tr><th>q<sub>gas,CH4</sub></th><th>q<sub>gas,H2</sub></th><th>q<sub>gas,CO</sub></th><th>q<sub>gas,CO2</sub></th><th>p<sub>gas,CH4</sub></th><th>p<sub>gas,H2</sub></th><th>p<sub>gas,CO</sub></th><th>p<sub>gas,CO2</sub></th></tr></thead><tbody><tr><td>R<sup>2</sup></td><td>0.95</td><td>0.90</td><td>0.90</td><td>NA</td><td>0.98</td><td>0.97</td><td>0.96</td><td>0.91</td></tr><tr><td>RMSE</td><td>0.29</td><td>0.38</td><td>0.27</td><td>NA</td><td>1.06</td><td>1.94</td><td>1.36</td><td>2.52</td></tr></tbody></table>

NA: Not available.

Fig. 5 displays the results of the simulations of the ADM1_ME with the estimated parameters against the experimental values. Regarding the experimental values, it is important to note that the behavior of the experimental data obtained from the two operational conditions is different. For OP1, there is an increase in the CH₄ content with a decrease in the CO₂ content which demonstrates the conversion of CO₂ into CH₄, however, for OP2, a decrease in both was observed. This decrease represents an inhibition of the biomethanation, which was due to the accumulation of acetate [54]. Even if different, the model fitted correctly both data sets.

<figure><img src="image_6.png" /></figure>
<caption>**Fig. 5.** Outlet gas flow rate ($q_{gas,i}$) and gas percent ($p_{gas,i}$) with the ADM1_ME. (A) OP1 and (B)

OP2.</caption>

For OP1 (Fig. 5-A), the model fitted correctly the experimental behavior of $q_{gas,CH4}$ and $q_{gas,CO}$.
However, it presents difficulties to reproduce the experimental behavior of $q_{gas,H2}$. This could be due
to unreported changes in the syngas composition, which here was assumed to be the same for the
simulations. Although, there is no experimental data available for $q_{gas,CO2}$, the model simulates
$q_{gas,CO2}$ whose dynamic is reducing progressively to compensate for the increase in methane content.
On the other hand, Fig. 5-B shows the ADM1_ME adjustment with the OP2, the model reproduces in
a better way the behavior of $p_{gas,CH4}$, $p_{gas,H2}$, $p_{gas,CO}$, and $p_{gas,CO2}$.

Table 7 reports the statistical evaluation of the ADM1_ME validation with the two operating conditions. Values higher than 0.74 and 0.82 were obtained in the $R^2$ with OP1 and OP2, except for $q_{gas,CH4}$. Concerning the RMSE, values less than 0.94 and 5.15 were exhibited with OP1 and OP2. The best predicted variables were $q_{gas,H2}$, $q_{gas,CO}$ ($R^2 > 0.74$ and RMSE < 0.94) with OP1. Concerning OP2, all the variables presented similar fitting ($R^2 > 0.82$ and RMSE < 5.15), resulting in a better model prediction.

<table><caption>**Table 6.** Statistical analysis for ADM1_ME validation with OP1 and OP2.</caption><thead><tr><th rowspan="2">Criteria</th><th colspan="4">OP1</th><th colspan="4">OP2</th></tr><tr><th>q<sub>gas,CH4</sub></th><th>q<sub>gas,H2</sub></th><th>q<sub>gas,CO</sub></th><th>q<sub>gas,CO2</sub></th><th>p<sub>gas,CH4</sub></th><th>p<sub>gas,H2</sub></th><th>p<sub>gas,CO</sub></th><th>p<sub>gas,CO2</sub></th></tr></thead><tbody><tr><td>R<sup>2</sup></td><td>0.39</td><td>0.74</td><td>0.81</td><td>NA</td><td>0.84</td><td>0.87</td><td>0.82</td><td>0.83</td></tr><tr><td>RMSE</td><td>0.94</td><td>0.94</td><td>0.22</td><td>NA</td><td>5.15</td><td>1.79</td><td>1.30</td><td>4.19</td></tr></tbody></table>

NA: Not available.

## 7.3. Model Analysis

Fig. 6 presents the concentration of $S_{liq,su}$, $S_{liq,bu}$, $S_{liq,pro}$, $S_{liq,ac}$, $S_{liq,CO}$, $S_{liq,H2}$, $S_{liq,CH4}$, and $S_{liq,CO2}$ with both operating conditions. The results were slightly similar in both cases, e.g., a final concentration of $5.9 \times 10^{-4}$ and $5.6 \times 10^{-4}$ gCOD/L were obtained for $S_{liq,su}$ with OP1 and OP2. Concentrations of $2.6 \times 10^{-2}$ and $2.5 \times 10^{-2}$ gCOD/L were obtained for $S_{liq,bu}$, while values of 0.15 and
0.16 gCOD/L were obtained for $S_{liq,pro}$. Steady state concentrations of 0.21 and 0.23 gCOD/L were obtained for $S_{liq,ac}$, 3.47×10⁻⁷ and 4.89×10⁻⁷ gCOD/L for $S_{liq,CO}$, whereas values of 1.74×10⁻⁷ and 1.92×10⁻⁷ gCOD/L were reached for $S_{liq,H₂}$. Concentrations of 6.84×10⁻² and 8.1×10⁻² gCOD/L were obtained for $S_{liq,CH₄}$. However, values of 1.11×10⁻⁵ and 3.20×10⁻³ mol/L were achieved for $S_{liq,CO₂}$ with the OP1 and OP2 2. This is evident from Fig. 5, in which CO₂ is almost exhausted for OP1 at stages IV-V but there is still some CO₂ for OP2 at stage IV.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Concentrations in the liquid phase: sugar ($S_{liq,su}$), butyrate ($S_{liq,bu}$), propionate ($S_{liq,pro}$), acetate ($S_{liq,ac}$), CO ($S_{liq,CO}$), H₂ ($S_{liq,H₂}$), CH₄ ($S_{liq,CH₄}$) and CO₂ ($S_{liq,CO₂}$). (A) OP1 and (B) OP2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Biomass concentrations that degrade the components: sugar ($X_{su}$), butyrate ($X_{bu}$), propionate ($X_{pro}$), acetate ($X_{ac}$), CO ($X_{CO}$), $H_2$ ($X_{H2}$). (A) OP1 and (B) OP2.</figcaption></figure>

Fig. 7 presents the biomass concentrations of $X_{su}$, $X_{bu}$, $X_{pro}$, $X_{ac}$, $X_{CO}$, and $X_{H2}$ which were simulated by the ADM1_ME with OP1 and OP2. The final biomass concentrations of $X_{su}$ were 0.62 and 3.56 gCOD/L for OP1 and OP2, respectively. Concentrations of 0.05 and 0.31 gCOD/L around were obtained for $X_{bu}$ and $X_{pro}$, whereas values of 0.25 and 1.18 gCOD/L were obtained for $X_{ac}$. At the end of the simulation, concentrations of 0.11 and 0.04 gCOD/L were obtained for $X_{CO}$, and 0.06 and 0.10 gCOD/L for $X_{H2}$ with OP1 and OP2. Regarding $X_{CO}$ and $X_{H2}$, it was observed that its concentration decreased drastically. However, it increased as the gas was added. This behavior illustrates the correlation between the added gas flow rate and the biomass concentration throughout the simulation of both OP1 and OP2. Results from OP1 indicate that the concentration of $X_{H2}$ fluctuated between 0.023 and 0.063 gCOD/L between stages II to V. Similarly, $X_{CO}$ values were found to range between 0.02 and 0.11 gCOD/L across the same stages. In contrast, OP2 displayed a
consistent $X_{H2}$ concentration of approximately 0.10 gCOD/L throughout all stages, while $X_{CO}$ values varied between 0.024 and 0.037 gCOD/L in stages I-V.

Fig. 8 displays the variables $S_{gas,H2}$, $S_{gas,CH4}$, $S_{gas,CO}$, and $S_{gas,CO2}$ simulated by the ADM1_ME with OP1 and OP2. Gas concentrations of 0.03 and 0.13 gCOD/L were obtained for $S_{gas,H2}$ with OP1 and OP2. The obtained gas concentrations at the end of the simulation for both operating conditions corresponded to 2.25 and 1.00 gCOD/L for $S_{gas,CH4}$, 0.03 and 0.10 gCOD/L for $S_{gas,CO}$, and 6.18×10⁻³ and 0.00 gCOD/L for $S_{gas,CO2}$.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Concentrations in the gas phase of H₂ ($S_{gas,H2}$), CH₄ ($S_{gas,CH4}$), CO ($S_{gas,CO}$) and CO₂ ($S_{gas,CO2}$). (A) OP1 and (B) OP2.</figcaption></figure>

The behavior for $S_{liq,j}$, $X_k$, and $S_{gas,i}$ between OP1 and OP2 are not directly comparable. They differ due to the dependence on either OLR or GLR. However, it is observed that the ADM1_ME reproduces different operational conditions and provides information about components such as $X_k$ that are difficult to measure.

Additionally, the model respects the COD balance which was calculated as proposed by Paudel et al. [61].

$$COD_{liq}^{in} + COD_{gas}^{in} + COD_{bio}^{in} = COD_{liq}^{out} + COD_{gas}^{out} + COD_{bio}^{out} \quad (36)$$
where $COD_{liq}^{in}$ and $COD_{gas}^{in}$ are the loading in the liquid and gas streams; $COD_{bio}^{in}$ is the initial COD biomass; $COD_{liq,out}$, $COD_{gas,out}$, and $COD_{bio,out}$ are the COD output in the liquid phase, the COD converted to produce biogas, and the assimilated COD for biomass growth, respectively.

For OP1, errors of 1.34, 2.88, 2.89, 2.45, 1.84, and 1.84% were obtained for the COD balance of the stages: reference, I, II, III, IV, and V, respectively. For OP2, the errors of the COD balance were 4.07, 3.90, 2.48, 2.07, and 1.14 for stages: reference, I, II, III, IV. The small errors could be due to the initial concentrations of biomass in the reactor.

Concerning the mass transfer, to improve the biomethanation process by increasing the syngas added, it is necessary to maintain the mass transfer capacity of the system. That is, in a BCR the system must be kept in a homogeneous bubbly flow regime, where the bubble size is nearly constant and dictated by the sparger design and the system properties [62]. In a CSTR, the impeller agitation speed must increase as the gas flow rates increases, ensuring high gas retention and complete dispersion (avoid flooding). In addition, in the case of aerated systems, the power consumption is lower than the non-aerated system due to the presence of cavities behind the agitator blades [63]. It is necessary to mention that this model is based on the hypothesis where both operational conditions are carried out satisfactorily: The BCR is operated in a homogeneous bubbly flow regime, and the power consumption in a CSTR is constant which allows maintaining the mass transfer.

As previously mentioned, one of the limiting factors of the performance of biomethanation is the mass transfer process [6,27,64]. The behavior is explored for both operating conditions OP1 and OP2. Fig. 9 displays the dynamic behavior of the volumetric mass transfer coefficients for OP1 and OP2. In both cases the values of $k_L a_i$ in the reference stage were 200 $1/d$, as proposed in [11]. However, $k_L a_i$ in OP1 (Fig. 9-A) depends only on the added gas flow rate. Therefore, their values decrease in stage I and then increased progressively up to 394.6, 347.1, 597.3, and 389.7 $1/d$ for $k_L a_{CO}$, $k_L a_{CH_4}$, $k_L a_{H_2}$, and $k_L a_{CO_2}$ in stage V. In contrast, the $k_L a_i$ in OP2 (Fig. 9-B) decreased in stage I below the value in the reference stage and increased progressively until reaching values of 56.22, 49.46, 85.08, and 55.54 $1/d$ for $k_L a_{CO}$, $k_L a_{CH_4}$, $k_L a_{H_2}$, and $k_L a_{CO_2}$, in stage V given the contribution provided by the gassed power.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Volumetric mass transfer coefficient of CH₄ (kLaCH₄), H₂ (kLaH₂), CO (kLaCO), and CO₂ (kLaCO₂). (A) OP1 and (B) OP2.</figcaption></figure>

## 8. Conclusions

An extension of the Anaerobic Digestion Model No. 1 (ADM1_ME) was proposed to represent the biomethanation process. The model extension was derived and assessed based on two operating conditions from the literature using two different bioreactor configurations: bubble column reactor and continuous stirred tank reactor, with two different substrates glucose and a mixture of primary sludge and ticked-disintegrated waste activated, respectively. Sensitivity analysis was performed by the Sobol' and Morris method to identify the candidate parameters to be estimated from which 14 out of 26 previously selected parameters (mainly kinetic parameters) presented high influence over the model outputs. Model validation was accurately performed where the variables $q_{gas,CO}$ and $q_{gas,H2}$ reported the best fitting ($R^2$ > 0.74 and RMSE < 0.94) for OP1, while for OP2 all the variables presented similar fitting ($R^2$ > 0.82 and RMSE < 5.15). Simulation results demonstrated that the proposed model could reproduce the gas outlet flow rates of the biomethanation process for BCR and CSTR while providing information about the dynamics of the biomasses involved in the process. Additionally, the model was able to simulate different operating conditions, and the use of various substrates, where an increase in CH₄ and a decrease in CO₂ content is expected. Further work will explore the application
of the model in model-based optimization to maximize yields and productivities of CH₄. Moreover, the model could be simulated for different conditions to generate data that could be used for machine-learning-based fault detection techniques.

## Credit authorship contribution statement

Juan C. Acosta-Pavas: Conceptualization, Investigation, Methodology, Software, Writing – original draft. Carlos. E Robles-Rodríguez: Conceptualization, Investigation, Methodology, Writing – review and editing. Jérôme Morchain: Conceptualization, Methodology, Writing – review and editing, Supervision, Project administration. Claire Dumas: Writing – review and editing. Arnaud Cockx: Writing – review and editing. César A. Aceves-Lara: Conceptualization, Investigation, Methodology, Writing – review and editing, Supervision, Project administration

## Acknowledgments

The authors would like to acknowledge the financial support of the Ministerio de Ciencias, Tecnología e Innovación (Minciencias) through Scholarship Program No. 860. This work has also benefited from a State grant managed by the National Research Agency under the "Investissements d'Avenir" program with the reference ANR-18-EURE-0021.

## References

[1] Brémond U, Bertrandias A, Steyer J-P, Bernet N, Carrere H. A vision of European biogas sector development towards 2030: Trends and challenges. Journal of Cleaner Production 2021;287:125065. https://doi.org/10.1016/j.jclepro.2020.125065.

[2] Gustafsson M, Svensson N, Eklund M, Fredriksson Möller B. Well-to-wheel climate performance of gas and electric vehicles in Europe. Transportation Research Part D: Transport and Environment 2021;97:102911. https://doi.org/10.1016/j.trd.2021.102911.

[3] Zupančič M, Možic V, Može M, Cimerman F, Golobič I. Current Status and Review of Waste-to-Biogas Conversion for Selected European Countries and Worldwide. Sustainability 2022;14:1823. https://doi.org/10.3390/su14031823.

[4] Iglesias R, Muñoz R, Polanco M, Díaz I, Susmozas A, Moreno AD, et al. Biogas from Anaerobic Digestion as an Energy Vector: Current Upgrading Development. Energies 2021;14:2742. https://doi.org/10.3390/en14102742.

[5] Rusmanis D, O'Shea R, Wall DM, Murphy JD. Biological hydrogen methanation systems – an overview of design and efficiency. Bioengineered 2019;10:604–34. https://doi.org/10.1080/21655979.2019.1684607.

[6] Rafrafi Y, Laguillaumie L, Dumas C. Biological Methanation of H2 and CO2 with Mixed Cultures: Current Advances, Hurdles and Challenges. Waste and Biomass Valorization 2020. https://doi.org/10.1007/s12649-020-01283-z.

[7] Grimalt-Alemany A, Asimakopoulos K, Skiadas IV, Gavala HN. Modeling of syngas biomethanation and catabolic route control in mesophilic and thermophilic mixed microbial consortia. Applied Energy 2020;262:114502. https://doi.org/10.1016/j.apenergy.2020.114502.

[8] Li C, Zhu X, Angelidaki I. Carbon monoxide conversion and syngas biomethanation mediated by different microbial consortia. Bioresource Technology 2020;314:123739. https://doi.org/10.1016/j.biortech.2020.123739.
703 [9] Mosey FE. Mathematical Modelling of the Anaerobic Digestion Process: Regulatory Mechanisms for the Formation of Short-Chain Volatile Acids from Glucose. Water Science and Technology 1983;15:209-32.
704 705 https://doi.org/10.2166/wst.1983.0168.

706 [10] Angelidaki I, Ellegaard L, Ahring BK. A comprehensive model of anaerobic bioconversion of complex substrates to biogas. Biotechnol Bioeng 1999;63:363-72. https://doi.org/10.1002/(SICI)1097-0290(19990505)63:3<363::AID-BIT13>3.0.CO;2-Z.

709 [11] Batstone DJ, Keller J, Angelidaki I, Kalyuzhnyi SV, Pavlostathis SG, Rozzi A, et al. The IWA Anaerobic Digestion Model No 1 (ADM1). Water Science and Technology 2002;45:65-73.
710 711 https://doi.org/10.2166/wst.2002.0292.

712 [12] Rosen C, Jeppsson U. Aspects on ADM1 Implementation within the BSM2 Framework. Department of Industrial Electrical Engineering and Automation (IEA), Lund University, Lund, Sweden: 2006.

714 [13] Bai J, Liu H, Yin B, Ma H, Chen X. Modified ADM1 for modeling free ammonia inhibition in anaerobic acidogenic fermentation with high-solid sludge. Journal of Environmental Sciences 2017;52:58-65.
715 716 https://doi.org/10.1016/j.jes.2016.03.004.

717 [14] Balde YM, Diop S, Tebbani S, Kante C. Modeling of a Continuous Anaerobic Digestion of Wastes. 2020 24th International Conference on System Theory, Control and Computing (ICSTCC), Sinaia, Romania: IEEE; 2020, p. 596-601. https://doi.org/10.1109/ICSTCC50638.2020.9259735.

720 [15] Blumensaat F, Keller J. Modelling of two-stage anaerobic digestion using the IWA Anaerobic Digestion Model No. 1 (ADM1). Water Research 2005;39:171-83. https://doi.org/10.1016/j.watres.2004.07.024.

722 [16] Sun H, Yang Z, Zhao Q, Kurbonova M, Zhang R, Liu G, et al. Modification and extension of anaerobic digestion model No.1 (ADM1) for syngas biomethanation simulation: From lab-scale to pilot-scale. Chemical Engineering Journal 2021;403:126177. https://doi.org/10.1016/j.cej.2020.126177.

725 [17] Ngu V, Morchain J, Cockx A. Spatio-temporal 1D gas-liquid model for biological methanation in lab scale and industrial bubble column. Chemical Engineering Science 2022;251:117478.
726 727 https://doi.org/10.1016/j.ces.2022.117478.

728 [18] Dar RA, Parmar M, Dar EA, Sani RK, Phutela UG. Biomethanation of agricultural residues: Potential, limitations and possible solutions. Renewable and Sustainable Energy Reviews 2021;135:110217.
729 730 https://doi.org/10.1016/j.rser.2020.110217.

731 [19] Laguillaumie L, Rafrafi Y, Moya-Leclair E, Delagnes D, Dubos S, Spérandio M, et al. Stability of ex situ biological methanation of H2/CO2 with a mixed microbial culture in a pilot scale bubble column reactor. Bioresource Technology 2022;354:127180. https://doi.org/10.1016/j.biortech.2022.127180.

734 [20] Czatzkowska M, Harnisz M, Korzeniewska E, Koniuszewska I. Inhibitors of the methane fermentation process with particular emphasis on the microbiological aspect: A review. Energy Science and Engineering 2020;8:1880-97. https://doi.org/10.1002/ese3.609.

737 [21] Chandra R, Takeuchi H, Hasegawa T. Methane production from lignocellulosic agricultural crop wastes: A review in context to second generation of biofuel production. Renewable and Sustainable Energy Reviews 2012;16:1462-76. https://doi.org/10.1016/j.rser.2011.11.035.

740 [22] Henze M, van Loosdrecht MCM, Ekama GA, Brdjanovic Damir. Biological Wastewater Treatment: Principles, Modeling and Design. IWA Publishing; 2019. https://doi.org/10.2166/9781780408613.

742 [23] Ashraf MT, Sieborg MU, Yde L, Rhee C, Shin SG, Triolo JM. Biomethanation in a thermophilic biotrickling filter — pH control and lessons from long-term operation. Bioresource Technology Reports 2020;11:100525. https://doi.org/10.1016/j.biteb.2020.100525.

745 [24] Bharathiraja B, Sudharsanaa T, Bharghavi A, Jayamuthunagai J, Praveenkumar R. Biohydrogen and Biogas - An overview on feedstocks and enhancement process. Fuel 2016;185:810-28.
746 747 https://doi.org/10.1016/j.fuel.2016.08.030.

748 [25] Saha S, Basak B, Hwang JH, Salama ES, Chatterjee PK, Jeon BH. Microbial Symbiosis: A Network towards Biomethanation. Trends in Microbiology 2020:1-17. https://doi.org/10.1016/j.tim.2020.03.012.

750 [26] Dev S, Saha S, Kurade MB, Salama ES, El-Dalatony MM, Ha GS, et al. Perspective on anaerobic digestion for biomethanation in cold environments. Renewable and Sustainable Energy Reviews 2019;103:85-95. https://doi.org/10.1016/j.rser.2018.12.034.

753 [27] Guiot SR, Cimpoia R, Carayon G. Potential of Wastewater-Treating Anaerobic Granules for Biomethanation of Synthesis Gas. Environ Sci Technol 2011;45:2006-12.
754 755 https://doi.org/10.1021/es102728m.

756 [28] Angelidaki I, Karakashev D, Batstone DJ, Plugge CM, Stams AJM. Biomethanation and its potential. vol. 494. 1st ed. Elsevier Inc.; 2011. https://doi.org/10.1016/B978-0-12-385112-3.00016-0.

758 [29] Liu R, Hao X, Wei J. Function of homoacetogenesis on the heterotrophic methane production with exogenous H2/CO2 involved. Chemical Engineering Journal 2016;284:1196-203.
759 760 https://doi.org/10.1016/j.cej.2015.09.081.
[30] Pan X, Zhao L, Li C, Angelidaki I, Lv N, Ning J, et al. Deep insights into the network of acetate metabolism in anaerobic digestion: focusing on syntrophic acetate oxidation and homoacetogenesis. Water Research 2021;190:116774. https://doi.org/10.1016/j.watres.2020.116774.

[31] Mulat DG, Mosbæk F, Ward AJ, Polag D, Greule M, Keppler F, et al. Exogenous addition of H2 for an in situ biogas upgrading through biological reduction of carbon dioxide into methane. Waste Management 2017;68:146-56. https://doi.org/10.1016/j.wasman.2017.05.054.

[32] Nauman EB. Chemical Reactor Design, Optimization, and Scaleup. Hoboken, NJ, USA: John Wiley & Sons, Inc.; 2008. https://doi.org/10.1002/9780470282076.

[33] Liu K, Phillips JR, Sun X, Mohammad S, Huhnke RL, Atiyeh HK. Investigation and modeling of gas-liquid mass transfer in a sparged and non-sparged continuous stirred tank reactor with potential application in syngas fermentation. Fermentation 2019;5. https://doi.org/10.3390/fermentation5030075.

[34] Gary, B. Tatterson. Fluid mixing and gas dispersion in agitated tanks 1991. https://doi.org/10.5860/choice.29-0927.

[35] Deckwer WD, Nguyen-Tien K, Kelkar BG, Shah YT. Applicability of axial dispersion model to analyze mass transfer measurements in bubble columns. AIChE Journal 1983;29:915-22. https://doi.org/10.1002/aic.690290607.

[36] Van't RK. Review of Measuring Methods and Results in Nonviscous Gas-Liquid Mass Transfer in Stirred Vessels 1979;18:357-64.

[37] Higbie R. The rate of absorption of a pure gas into still liquid during short periods of exposure. vol. 31. Transaction of the American Institute of Chemical Engineers (AIChe); 1935.

[38] Jensen MB, Ottosen LDM, Kofoed MVW. H2 gas-liquid mass transfer: A key element in biological Power-to-Gas methanation. Renewable and Sustainable Energy Reviews 2021;147:111209. https://doi.org/10.1016/j.rser.2021.111209.

[39] Damblin G, Couplet M, Iooss B. Numerical studies of space-filling designs: optimization of Latin Hypercube Samples and subprojection properties. Journal of Simulation 2013;7:276-89. https://doi.org/10.1057/jos.2013.16.

[40] Sepulveda FD, Cisternas LA, Gálvez ED. Global sensitivity analysis of a mineral processing flowsheet. Computer Aided Chemical Engineering, vol. 32, Elsevier; 2013, p. 913-8. https://doi.org/10.1016/B978-0-444-63234-0.50153-6.

[41] Sohier H, Farges J-L, Piet-Lahanier H. Improvement of the Representativity of the Morris Method for Air-Launch-to-Orbit Separation. IFAC Proceedings Volumes 2014;47:7954-9. https://doi.org/10.3182/20140824-6-ZA-1003.01968.

[42] Tosin M, Côrtes AMA, Cunha A. A Tutorial on Sobol' Global Sensitivity Analysis Applied to Biological Models. In: da Silva FAB, Carels N, Trindade dos Santos M, Lopes FJP, editors. Networks in Systems Biology, vol. 32, Cham: Springer International Publishing; 2020, p. 93-118. https://doi.org/10.1007/978-3-030-51862-2_6.

[43] Feng K, Lu Z, Yang C. Enhanced Morris method for global sensitivity analysis: good proxy of Sobol' index. Struct Multidisc Optim 2019;59:373-87. https://doi.org/10.1007/s00158-018-2071-7.

[44] Kucherenko S, Delpuech B, Iooss B, Tarantola S. Application of the control variate technique to estimation of total sensitivity indices. Reliability Engineering & System Safety 2015;134:251-9. https://doi.org/10.1016/j.ress.2014.07.008.

[45] Zhang X, Trame M, Lesko L, Schmidt S. Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models. CPT Pharmacometrics Syst Pharmacol 2015;4:69-79. https://doi.org/10.1002/psp4.6.

[46] Morio J. Global and local sensitivity analysis methods for a physical system. Eur J Phys 2011;32:1577-83. https://doi.org/10.1088/0143-0807/32/6/011.

[47] Ochoa MP, Estrada V, Di Maggio J, Hoch PM. Dynamic global sensitivity analysis in bioreactor networks for bioethanol production. Bioresource Technology 2016;200:666-79. https://doi.org/10.1016/j.biortech.2015.10.069.

[48] Iooss B, Lemaître P. A Review on Global Sensitivity Analysis Methods. In: Dellino G, Meloni C, editors. Uncertainty Management in Simulation-Optimization of Complex Systems, vol. 59, Boston, MA: Springer US; 2015, p. 101-22. https://doi.org/10.1007/978-1-4899-7547-8_5.

[49] Sepúlveda FD, Cisternas LA, Gálvez ED. The use of global sensitivity analysis for improving processes: Applications to mineral processing. Computers & Chemical Engineering 2014;66:221-32. https://doi.org/10.1016/j.compchemeng.2014.01.008.

[50] Sobol' IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation 2001;55:271-80. https://doi.org/10.1016/S0378-4754(00)00270-6.
819 [51] Morris MD. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics 1991;33:161-74. https://doi.org/10.1080/00401706.1991.10484804.

820 [52] Saltelli A, editor. Sensitivity analysis in practice: a guide to assessing scientific models. Reprinted.
822 Chichester Weinheim: Wiley; 2007.

823 [53] Saltelli A, editor. Global sensitivity analysis: the primer. Chichester, West Sussex: Wiley; 2008.

824 [54] Andreides D, Pokorna D, Zabranska J. Assessing the syngas biomethanation in anaerobic sludge digestion under different syngas loading rates and homogenisation. Fuel 2022;320:123929.
825 https://doi.org/10.1016/j.fuel.2022.123929.
826

827 [55] Ahnert M, Schalk T, Brückner H, Effenberger J, Kuehn V, Krebs P. Organic matter parameters in WWTP – a critical review and recommendations for application in activated sludge modelling. Water Science and Technology 2021;84:2093–112. https://doi.org/10.2166/wst.2021.419.

830 [56] Batstone DJ, Balthes C, Barr K. Model assisted startup of anaerobic digesters fed with thermally
831 hydrolysed activated sludge. Water Science 2010:6.

832 [57] Velez S. Carlos M. Global sensitivity and uncertainty analysis (GSUA)(https://www.mathworks.com/matlabcentral/fileexchange/47758-global-sensitivity-and-uncertainty-analysis-gsua), 2022.
833
834

835 [58] Mr. Sensitivity Analysis - Morris method (advanced) (https://www.mathworks.com/matlabcentral/fileexchange/48884-sensitivity-analysis-morris-method-advanced), 2022.
836
837

838 [59] Rodriguez-Fernandez M, Kucherenko S, Pantelides C, Shah N. Optimal experimental design based on global sensitivity analysis. Computer Aided Chemical Engineering, vol. 24, Elsevier; 2007, p. 63–8.
839
840 https://doi.org/10.1016/S1570-7946(07)80034-4.

841 [60] Asprey SP, Macchietto S. Statistical tools for optimal dynamic model building. Computers & Chemical Engineering 2000;24:1261-7. https://doi.org/10.1016/S0098-1354(00)00328-8.
842

843 [61] Paudel S, Kang Y, Yoo Y-S, Seo GT. Hydrogen Production in the Anaerobic Treatment of Domestic-
844 Grade Synthetic Wastewater. Sustainability 2015;7:16260-72. https://doi.org/10.3390/su71215814.

845 [62] Hissanaga AM, Padoin N, Paladino EE. Mass transfer modeling and simulation of a transient homogeneous bubbly flow in a bubble column. Chemical Engineering Science 2020;218:115531. https://doi.org/10.1016/j.ces.2020.115531.

848 [63] Gabelle J -C., Augier F, Carvalho A, Rousset R, Morchain J. Effect of tank size on kLa and mixing time in aerated stirred reactors with non-newtonian fluids. Can J Chem Eng 2011;89:1139-53. https://doi.org/10.1002/cjce.20571.

851 [64] Luo G, Wang W, Angelidaki I. Anaerobic Digestion for Simultaneous Sewage Sludge Treatment and CO
852 Biomethanation: Process Performance and Microbial Ecology. Environ Sci Technol
853 2013:130904143045005. https://doi.org/10.1021/es401018d.

854