Mathematical Model-Based Optimization of Trace Metal Dosage in Anaerobic Batch Bioreactors

Abstract

Anaerobic digestion (AD) is a promising and yet a complex waste-to-energy technology. To optimize such a process, precise modeling is essential. Developing complex, mechanistically inspired AD models can result in an overwhelming number of parameters that require calibration. This study presents a novel approach that considers the role of trace metals (Ca, K, Mg, Na, Co, Cr, Cu, Fe, Ni, Pb, and Zn) in the modeling, numerical simulation, and optimization of the AD process in a batch bioreactor. In this context, BioModel is enhanced by incorporating the influence of metal activities on chemical, biochemical, and physicochemical processes. Trace metal-related parameters are also included in the calibration of all model parameters. The model’s reliability is rigorously validated by comparing simulation results with experimental data. The study reveals that perturbations of 5% in model parameter values significantly increase the discrepancy between simulated and experimental results up to threefold. Additionally, the study highlights how precise optimization of metal additives can enhance both the quantity and quality of biogas production. The optimal concentrations of trace metals increased biogas and CH₄ production by 5.4% and 13.5%, respectively, while H₂, H₂S, and NH₃ decreased by 28.2%, 43.6%, and 42.5%, respectively.

1. Introduction

Ever increasing demands for reducing waste and greenhouse gas emissions coupled with increasing renewable energy sources are the driving force for the development of circular waste management technologies such as anaerobic digestion (AD) technology. AD enables reduction of waste entering landfills and has the potential to reduce greenhouse gas (GHG) emissions [1,2,3,4]; the World Biogas Association confirmed in a recent report that AD technology has a potential to reduce GHG emissions in the range between 3290 and 4360 Mt CO_2,eq, which is equivalent to 10–13% of the world’s current emissions [4]. Furthermore, the production of biogas by AD processes has reached approximately 66 billion m³ in 2020, which is equivalent to an energy content of 1.52 EJ. Interestingly, the International Renewable Energy Agency (IRENA) forecasted a tenfold increase in biogas pathways, leading to a supply of 15 EJ by 2050 as part of the Transforming Energy Scenario [5].

The energy, environmental, and financial incentives associated with setting up an AD system necessitates rigorous engineering analysis of the process, where microorganisms in a bioreactor breakdown organic matter (OM) in the absence of oxygen to biogas. The quality and quantity of biogas production depend mainly on AD process conditions, feeding strategy, and feedstock substrate composition (macro and micronutrients). In addition to the microbial community, enzymes and organic degradable matters, such as carbohydrates, proteins, lipids, inorganic compounds, and trace metals (TMs), play a very important role. Although the AD system has been extensively studied, the role of TMs is still under scrutiny and has not sufficiently undergone rigorous engineering analysis to establish tradeoffs between efficiency and cost.

Light metals (Ca, K, Mg, Na, etc.) and heavy metals (Co, Cr, Cu, Fe, Ni, Pb, Zn, etc.) are integral to the biochemical processes occurring in the bioreactor [6,7,8] that facilitate biodegradation. During biodegradation, TMs influence the biogas production quality and quantity by impacting different stages: hydrolysis, acidogenesis, acetogenesis, and methanogenesis. The impact of the TMs on these processes is a consequence of the effect of TMs on microbial consortium activity. Each individual microorganism is affected by various TMs in multiple ways. For example, Co supports the growth and activity of acetogenic bacteria, while Fe, Zn, and Ni are important methanogens building elements [9]. The involvement of TMs (Fe, Ni, Cu, Ca, Mg, and Co) into chemical and physicochemical processes, where sulfide, carbonate, and phosphate precipitates are formed, influences the pH value during the AD process and also the quality and quantity of the produced biogas. Fe can reduce H₂S in the produced biogas by shifting the chemical equilibrium toward the formation of sulfide precipitates [10,11,12]. Likewise, a Zn concentration below 470 mg L⁻¹ can enhance the conversion of organic matter to CH₄ because zinc ions are a critical element of the active sites in multiple enzymes [13,14]. In general, a lack of TMs reduces microbial growth and activity, which significantly limits the AD process, resulting in accumulation of metabolic intermediates [15]. One intermediate is volatile fatty acids, which cause bioreactor acidification and consequently reduce methane yield. Since some substrates for the AD process may lack adequate amounts of TMs, TM supplementation is necessary in order to maintain an efficient AD process.

Contrarily, many investigations have concluded that an excess of essential TMs may also strongly inhibit the AD process [14,15,16]. It is worth noting that the reported inhibitory concentration of each metal in a bioreactor depends on the interaction of various AD parameters, especially the available microbial community, AD process conditions, and feedstock substrate composition [17]. For example, Chen et al. [18] reported that a concentration of Na over 3.5 g L⁻¹ is inhibitory to methanogens at mesophilic temperature, while Geng et al. [19] reported that a Na concentration over 0.35 g L⁻¹ can cause hyperosmotic stress, which results in the cellular dehydration of methanogens and a decrease in dehydrogenase activity during food waste AD. Abdel-Shafy and Mansour [20] reported that a Cr concentration over 0.775 g L⁻¹ has an inhibitory effect on AD of sludge, while Jha and Schmidt [21] reported that a Cr concentration of 0.012 g L⁻¹ may reduce methanogens by 50%. Zhang et al. [16] reported that Fe concentrations over 5.65 g L⁻¹ can have inhibitory effects on the AD process of cow dung and Phragmites straw. On the other hand, Wang et al. [22] showed that a concentration of Fe above 2.9 g L⁻¹ in the AD of waste-activated sludge can decrease CH₄ production due to the formation of carbonate precipitates, thereby reducing the ${\mathrm{C}\mathrm{O}}_3^{2-}$ and ${\mathrm{H}\mathrm{C}\mathrm{O}}_3^{-}$ content in the bulk phase and decreasing the availability of CO₂ for hydrogen-utilizing methanogens. Furthermore, Xing et al. [14] investigated the long- and short-term effects of the Zn concentration during excess activated sludge AD and concluded that a Zn concentration over 670 mg L⁻¹ reduced methanogenic activities markedly, indicating that an excessive concentration of Zn significantly impairs system stability. These results show that optimal concentrations of metals can vary for various AD scenarios.

As the experiments related to AD are time-consuming, labor intensive, and expensive, numerical simulations (NSs) provide a cost-effective solution when there is a large number of operating parameters with complex interrelationships. So far, NS, which provides a useful tool for AD process understanding and optimization, has been performed using more or less complex mechanistically inspired mathematical models [6,15,23,24,25,26,27,28,29,30,31,32,33,34,35]. It is clear that with the complexity of AD models, which consider various biochemical, chemical, and physicochemical processes, the number of model parameters increases. Since these parameters are unknown or hard to determine, various efficient procedures for model parameters calibration have been developed in multiple studies [28,29,30,35]. However, many researchers calibrate only the most important model parameter, while the estimated values from the literature are typically used for the less important parameters. This simplification neglects the interactions among the less significant model parameters [28,29,30]. Some of these model parameters are related to the activities of TMs, which were included to some extent in some complex AD models [15,26,28,29,30,32,33,34]. For example, Frunzo et al. [15] and Flores-Alsina et al. [26] included the interactions of iron with phosphorus and sulfur into the modeling of the AD process in a CSTR bioreactor. For NS of the AD process in a batch bioreactor, Maharaj et al. [32,33,34] upgraded the AD model by including the activities of Ca, Co, Fe, Mg, and Ni in physicochemical processes; all TM-related model parameters were estimated according to the literature. Kegl et al. [28,29,30] upgraded BioModel for NS of AD performance in a CSTR bioreactor by considering TMs activities; all TM-related parameters were calibrated. However, challenges remain in calibrating and validating such complex models due to the lack of reported data on the range of TM concentrations required for optimal microbial activity, particularly in relation to methane production and other key intermediates in the AD process [34]. Therefore, there is still a need for adequate modeling of the role and activities of all key TMs during the AD process and adequate calibration of all model parameters, including those related to TMs. To the author’s knowledge, no study has yet addressed the importance of accurately calibrating all model parameters. Additionally, there is still a lack of sufficient numerical optimization regarding the optimal concentrations of TMs in the influent feedstock.

This paper investigates the impact of TMs on AD performance, focusing on the modeling, calibration, validation, and numerical optimization of TM concentrations in lab-scale batch bioreactors, with the aim of improving both the understanding and practical application of TM management in AD processes. To explore the impact of various TMs on the AD process, a self-developed BioModel is enhanced to simulate the dynamics of metal interactions throughout the process. The enhanced BioModel incorporates 187 model parameters, all of which require careful calibration. To assess the precision of this calibration, both low and high perturbations of the parameters are thoroughly investigated and analyzed. The calibrated and validated BioModel, tested under diverse operating conditions across different bioreactors, is then integrated into an optimization procedure using a gradient-based approximation method. The optimization aims to determine the optimal quantities of metals added or removed in the influent feedstock, defined using a range of multi-objective and constraint functions. Finally, the results of the BioModel calibration, validation, and numerical optimization are presented and discussed in detail.

2. Materials and Methods

At first, the experimental data of the AD process are presented, followed by a detailed description of the enhanced BioModel. In the following, the paper outlines the procedures employed for BioModel calibration and the optimization process, which aims to fine-tune the quantities of trace metals (TMs) in the feedstock to enhance the AD process, using a sophisticated gradient-based approximation method.

2.1. Experimental Data

The AD processes of cattle and cow manure in three-batch-mode lab-scale bioreactors (B1, B2, B3) were observed. The bioreactor data are given in

Table 1. Bioreactor data.

	Bioreactor
Parameter	B1	B2	B3
$\mathrm{pH} \mathrm{value}, \mathrm{p}\mathrm{H} (/)$	8.00	8.00	7.4
$\mathrm{Temperature}, T (℃)$	35.00	37.00	38.5
$\mathrm{Liquid} \mathrm{volume} \mathrm{of} \mathrm{bioreactor}, V_{\mathrm{l}\mathrm{i}\mathrm{q}}\left(\mathrm{L}\right)$	1.00	0.25	4.5
$\mathrm{Gas} \mathrm{volume} \mathrm{of} \mathrm{bioreactor}, V_{\mathrm{g}}\left(\mathrm{L}\right)$	0.10	0.20	3.5
$\mathrm{Total} \mathrm{pressure} \mathrm{in} \mathrm{bioreactor}, p_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}\left(\mathrm{b}\mathrm{a}\mathrm{r}\right)$	1.0059	1.0059	1.0073

. The animal manure is composed of carbohydrates, proteins, lipids, inert matters, inorganic carbon, inorganic nitrogen, inorganic phosphorus, inorganic potassium, inorganic sulfur in the form of sulfide and sulfate, and other elements and compounds [27,36,37,38]. The input data of the feedstock composition for the AD process in bioreactors B1 [27,36], B2 [37], and B3 [38] are given in

Table 2. Input data of feedstock composition.

Parameter	Value
$\mathrm{Carbohydrate} \mathrm{concentration}, c_{\mathrm{c}\mathrm{h}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	15.0000
$\mathrm{Protein} \mathrm{concentration}, c_{\mathrm{p}\mathrm{r}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	35.0000
$\mathrm{Lipid} \mathrm{concentration}, c_{\mathrm{l}\mathrm{i}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	3.0000
$\mathrm{Inert} \mathrm{matter} \mathrm{concentration}, c_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{t}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	50.0000
$\mathrm{Inorganic} \mathrm{carbon} \mathrm{concentration}, c_{\mathrm{C}\mathrm{i}\mathrm{o}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	2.5000
$\mathrm{Inorganic} \mathrm{nitrogen} \mathrm{concentration}, c_{\mathrm{N}\mathrm{i}\mathrm{o}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	2.0000
$\mathrm{Inorganic} \mathrm{phosphorous} \mathrm{concentration}, c_{\mathrm{P}\mathrm{i}\mathrm{o}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	2.6000
$\mathrm{Inorganic} \mathrm{sulfur} \left(\mathrm{sulfide}\right) \mathrm{concentration}, c_{\mathrm{S}\mathrm{i}\mathrm{o}\mathrm{i}\mathrm{d}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.1000
$\mathrm{Inorganic} \mathrm{sulfur} \left(\mathrm{sulfate}\right) \mathrm{concentration}, c_{\mathrm{S}\mathrm{i}\mathrm{o}\mathrm{a}\mathrm{t}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.5000
$\mathrm{Nitrogen} \mathrm{dioxide} \mathrm{concentration}, c_{{\mathrm{N}\mathrm{O}}_2} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0080
$\mathrm{Calcium} \mathrm{concentration}, c_{\mathrm{C}\mathrm{a}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	3.0000
$\mathrm{Cobalt} \mathrm{concentration}, c_{\mathrm{C}\mathrm{o}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0200
$\mathrm{Chromium} \mathrm{concentration}, c_{\mathrm{C}\mathrm{r}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0030
$\mathrm{Copper} \mathrm{concentration}, c_{\mathrm{C}\mathrm{u}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0100
$\mathrm{Iron} \mathrm{concentration}, c_{\mathrm{F}\mathrm{e}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.1100
$\mathrm{Potassium} \mathrm{concentration}, c_{\mathrm{K}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	3.2000
$\mathrm{Magnesium} \mathrm{concentration}, c_{\mathrm{M}\mathrm{g}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.8000
$\mathrm{Sodium} \mathrm{concentration}, c_{\mathrm{N}\mathrm{a}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.3000
$\mathrm{Nickel} \mathrm{concentration}, c_{\mathrm{N}\mathrm{i}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0010
$\mathrm{Lead} \mathrm{concentration}, c_{\mathrm{P}\mathrm{b}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0002
$\mathrm{Zinc} \mathrm{concentration}, c_{\mathrm{Z}\mathrm{n}} \left({\mathrm{g} \mathrm{L}}^{-1}\right)$	0.0500

.

2.2. Numerical Simulation of the AD Process

To ensure a reliable and efficient NS of AD performance in a batch-mode AD process, an upgraded BioModel was developed and incorporated into the active set optimization (ASO) procedure.

Enhanced BioModel. The upgraded BioModel is based on the BioModel presented in [28,29,30], which considers chemical (acid–base reactions), physicochemical (mass transfer from liquid to gas phase and the precipitation), and biochemical (hydrolysis, acidogenesis, acetogenesis, and methanogenesis) processes. In the biochemical reactions, the inhibitory effects of various compounds and TMs on bacteria growth are modeled with non-competitive inhibition.

The main modifications are related to the modeling of the activities of various TMs from the set $I_{\mathrm{T}\mathrm{M}\mathrm{s}}=\left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{C}\mathrm{l},\mathrm{C}\mathrm{r},\mathrm{F}\mathrm{e},\mathrm{K},\mathrm{M}\mathrm{g},\mathrm{N}\mathrm{a},\mathrm{N}\mathrm{i},\mathrm{P}\mathrm{b},\mathrm{S}\mathrm{e},\mathrm{Z}\mathrm{n}\right\}$ as these TMs are present in the feeding substrate. The effects of all TMs during the AD process are incorporated into the upgraded BioModel using Equations (1)–(8).

At first, the effects of TMs on the pH value in the bioreactor are determined using Equation (1). The ${\mathrm{H}}^{+}$ concentration derivative, needed for pH determination, can be expressed as follows:

$${\displaystyle \frac{d\left[{\mathrm{H}}^{+}\right]}{dt}}={\displaystyle \frac{d\left[{\mathrm{H}}^{+}\right]}{d{C}_{\mathrm{h}}}} {\displaystyle \frac{d{C}_{\mathrm{h}}}{dt}}$$

(1)

Here, $C_{\mathrm{h}}$ denotes the sum of ionic metals concentrations obtained from the following charge balance equation:

$$\begin{array}{ll}{C}_{\mathrm{h}}=\left[{\mathrm{p}\mathrm{r}\mathrm{o}}^{-}\right]& +\left[{\mathrm{b}\mathrm{u}}^{-}\right]+\left[{\mathrm{v}\mathrm{a}}^{-}\right]+\left[{\mathrm{a}\mathrm{c}}^{-}\right]-\left[{\mathrm{N}\mathrm{H}}_4^{+}\right]+\left[\mathrm{H}{\mathrm{C}\mathrm{O}}_3^{-}\right]+2\left[{\mathrm{C}\mathrm{O}}_3^{2-}\right]+\left[\mathrm{H}{\mathrm{S}\mathrm{O}}_4^{-}\right]\\ & +2\left[{\mathrm{S}\mathrm{O}}_4^{2-}\right]+\left[{\mathrm{H}\mathrm{S}}^{-}\right]+2\left[{\mathrm{S}}^{2-}\right]+\left[{\mathrm{H}}_2{\mathrm{P}\mathrm{O}}_4^{-}\right]+2\left[{\mathrm{H}\mathrm{P}\mathrm{O}}_4^{2-}\right]\\ & +3\left[{\mathrm{P}\mathrm{O}}_4^{3-}\right]+\left[{\mathrm{N}\mathrm{O}}_2^{-}\right]+\left[{\mathrm{A}\mathrm{n}}^{-}\right]-2\left[{\mathrm{C}\mathrm{a}}^{2+}\right]-\left[{\mathrm{C}\mathrm{a}\mathrm{t}}^{+}\right]+\left[{\mathrm{C}\mathrm{l}}^{-}\right]\\ & -2\left[{\mathrm{C}\mathrm{o}}^{2+}\right]-2\left[{\mathrm{C}\mathrm{r}}^{2+}\right]-2\left[{\mathrm{C}\mathrm{u}}^{2+}\right]-2\left[{\mathrm{F}\mathrm{e}}^{2+}\right]-\left[{\mathrm{K}}^{+}\right]-2\left[{\mathrm{M}\mathrm{g}}^{2+}\right]\\ & -\left[{\mathrm{N}\mathrm{a}}^{+}\right]-2\left[{\mathrm{N}\mathrm{i}}^{2+}\right]-2\left[{\mathrm{P}\mathrm{b}}^{2+}\right]-2\left[{\mathrm{S}\mathrm{e}}^{2+}\right]-2\left[{\mathrm{Z}\mathrm{n}}^{2+}\right]\end{array}$$

The role of TMs in precipitation is included into the modeling of liquid–solid precipitation rates and the concentrations of precipitates from the set $I_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}}=\{{\mathrm{C}\mathrm{a}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{C}\mathrm{o}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{C}\mathrm{u}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{F}\mathrm{e}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{M}\mathrm{g}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{N}\mathrm{i}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{P}\mathrm{b}\mathrm{C}\mathrm{O}}_3$, ${\mathrm{Z}\mathrm{n}\mathrm{C}\mathrm{O}}_3$, $\mathrm{F}\mathrm{e}\mathrm{S}, \mathrm{C}\mathrm{o}\mathrm{S}$, $\mathrm{C}\mathrm{u}\mathrm{S}, \mathrm{N}\mathrm{i}\mathrm{S}$, $\mathrm{P}\mathrm{b}\mathrm{S}, \mathrm{Z}\mathrm{n}\mathrm{S}$, ${\mathrm{C}\mathrm{a}}_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2$, ${\mathrm{C}\mathrm{o}}_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2$, ${\mathrm{F}\mathrm{e}}_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2$, ${\mathrm{N}\mathrm{i}}_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2$, $\mathrm{M}\mathrm{g}{{\mathrm{N}\mathrm{H}}_4\mathrm{P}\mathrm{O}}_4$, $\mathrm{K}\mathrm{M}\mathrm{g}{\mathrm{P}\mathrm{O}}_4\}$. The precipitation rates of sulfides ${\rho }_{{\mathrm{l}-\mathrm{s},\mathrm{S}}^{2-}} \left(\mathrm{g} {\mathrm{L}}^{-1}{\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right),$ carbonates ${\rho }_{\mathrm{l}-\mathrm{s}{,\mathrm{C}\mathrm{O}}_3^{2-}} \left(\mathrm{g} {\mathrm{L}}^{-1}{\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right),$ phosphates ${\rho }_{\mathrm{l}{-\mathrm{s},\mathrm{P}\mathrm{O}}_4^{3-}} \left(\mathrm{g} {\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$, and ionic forms of TMs ${\rho }_{\mathrm{l}-\mathrm{s},i} \left(\mathrm{g} {\mathrm{L}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right),i\in I_{\mathrm{T}\mathrm{M}\mathrm{s}}$ are determined as follows:

$${\rho }_{{\mathrm{l}-\mathrm{s},\mathrm{S}}^{2-}}=\sum _i\begin{array}{c}{k}_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i\mathrm{S}}{\left({\left(c_{i^{2+}}{c}_{{\mathrm{S}}^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(K_{\mathrm{s}\mathrm{p},i\mathrm{S}}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2\\ \end{array}, i\in \left\{\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{N}\mathrm{i},\mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}\right\}$$

(2)

$${\rho }_{\mathrm{l}-\mathrm{s}{,\mathrm{C}\mathrm{O}}_3^{2-}}=\sum _i{k}_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i{\mathrm{C}\mathrm{O}}_3}{\left({\left(c_{i^{2+}}{c}_{{\mathrm{C}\mathrm{O}}_3^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(K_{\mathrm{s}\mathrm{p},i{\mathrm{C}\mathrm{O}}_3}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{M}\mathrm{g},\mathrm{N}\mathrm{i}, \mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}\right\}$$

(3)

$$\begin{array}{ll}{\rho }_{\mathrm{l}{-\mathrm{s},\mathrm{P}\mathrm{O}}_4^{3-}}& =2{\displaystyle \sum _i}{k}_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i_2{\left({\mathrm{P}\mathrm{O}}_4\right)}_3}{\left({\left({\left(c_{i^{2+}}\right)}^3{\left(c_{{\mathrm{P}\mathrm{O}}_4^{3-}}\right)}^2\right)}^{{\displaystyle \frac{1}{5}}}-{\left(K_{\mathrm{s}\mathrm{p},i_2{\left({\mathrm{P}\mathrm{O}}_4\right)}_3}\right)}^{{\displaystyle \frac{1}{5}}}\right)}^2\\ & +k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{M}\mathrm{g}{\mathrm{N}\mathrm{H}}_4{\mathrm{P}\mathrm{O}}_4}{\left({\left(c_{{\mathrm{M}\mathrm{g}}^{2+}}{c}_{{\mathrm{N}\mathrm{H}}_4^{+}}{c}_{{\mathrm{P}\mathrm{O}}_4^{3-}}\right)}^{{\displaystyle \frac{1}{3}}}-{\left(K_{\mathrm{s}\mathrm{p},\mathrm{M}\mathrm{g}{\mathrm{N}\mathrm{H}}_4{\mathrm{P}\mathrm{O}}_4}\right)}^{{\displaystyle \frac{1}{3}}}\right)}^2\\ & +k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{K}\mathrm{M}\mathrm{g}{\mathrm{P}\mathrm{O}}_4}{\left({\left(c_{{\mathrm{K}}^{+}}{c}_{{\mathrm{M}\mathrm{g}}^{2+}}{c}_{{\mathrm{P}\mathrm{O}}_4^{3-}}\right)}^{{\displaystyle \frac{1}{3}}}-{\left(K_{\mathrm{s}\mathrm{p},\mathrm{K}\mathrm{M}\mathrm{g}{\mathrm{P}\mathrm{O}}_4}\right)}^{{\displaystyle \frac{1}{3}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{F}\mathrm{e},\mathrm{N}\mathrm{i}\right\}\end{array}$$

(4)

$$\begin{array}{ll}{\rho }_{\mathrm{l}-\mathrm{s},i^{2+}}& =k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},{i\mathrm{C}\mathrm{O}}_3} {\left({\left(c_{i^{2+}} c_{{\mathrm{C}\mathrm{O}}_3^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(K_{\mathrm{s}\mathrm{p},{i\mathrm{C}\mathrm{O}}_3}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2+{3 k}_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2}{\left({\left({\left(c_{i^{2+}}\right)}^3 {\left(c_{{\mathrm{P}\mathrm{O}}_4^{3-}}\right)}^2\right)}^{{\displaystyle \frac{1}{5}}}-{\left(K_{\mathrm{s}\mathrm{p},i_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2}\right)}^{{\displaystyle \frac{1}{5}}}\right)}^2\\ & +k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i\mathrm{S}} {\left({\left(c_{i^{2+}} c_{{\mathrm{S}}^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(K_{\mathrm{s}\mathrm{p},i\mathrm{S}}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{K},\mathrm{M}\mathrm{g},\mathrm{N}\mathrm{i},\mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}\right\}\end{array}$$

(5)

where the symbols $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t}} \left({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$ and $K_{\mathrm{s}\mathrm{p}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ are the crystallization rate and solubility product constants of precipitates from the set $I_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}}$, respectively, while the symbol $c_k \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ denotes the concentration of the corresponding $k^{\mathrm{t}\mathrm{h}}$ compound in ionic form (sulfide, carbonate, phosphate, and TMs).

The precipitation rates of sulfides in Equation (2) are obtained from precipitation/dissolution reactions ${\left(\mathrm{T}\mathrm{M}\mathrm{s}\right)}_{\left(\mathrm{l}\right)}^{2+}+{\mathrm{S}}_{\left(\mathrm{l}\right)}^{2-}\leftrightarrow {\left(\mathrm{T}\mathrm{M}\mathrm{s}\right)\mathrm{S}}_{\left(\mathrm{s}\right)}$, where the TMs are Co, Cu, Fe, Ni, Pb, and Zn. The precipitation/dissolution reactions ${\left(\mathrm{T}\mathrm{M}\mathrm{s}\right)}_{\left(\mathrm{l}\right)}^{2+}+{\mathrm{C}\mathrm{O}}_{3 \left(\mathrm{l}\right)}^{2-}\leftrightarrow {\left(\mathrm{T}\mathrm{M}\mathrm{s}\right)\mathrm{C}\mathrm{O}}_{3 \left(\mathrm{s}\right)}$, where the TMs are Ca, Co, Cu, Fe, Mg, Ni, Pb, and Zn, are the basis for precipitation rates of carbonates in Equation (3). Furthermore, Equation (4) is obtained from the reaction ${3 \left(\mathrm{T}\mathrm{M}\mathrm{s}\right)}_{\left(\mathrm{l}\right)}^{2+}+{2\mathrm{P}\mathrm{O}}_{4 \left(\mathrm{l}\right)}^{3-}\leftrightarrow {{\left(\mathrm{T}\mathrm{M}\mathrm{s}\right)}_3\left({\mathrm{P}\mathrm{O}}_4\right)}_{2 \left(\mathrm{s}\right)}$, where the TMs are Ca, Co, Fe, and Ni, and from the reactions ${\mathrm{M}\mathrm{g}}_{\left(\mathrm{l}\right)}^{2+}+{\mathrm{N}\mathrm{H}}_{4 \left(\mathrm{l}\right)}^{+}{+\mathrm{P}\mathrm{O}}_{4 \left(\mathrm{l}\right)}^{3-}\leftrightarrow {\mathrm{M}\mathrm{g}{\mathrm{N}\mathrm{H}}_4\mathrm{P}\mathrm{O}}_{4 \left(\mathrm{s}\right)}$ and ${\mathrm{K}}_{\left(\mathrm{l}\right)}^{+}+{\mathrm{M}\mathrm{g}}_{\left(\mathrm{l}\right)}^{2+}{+\mathrm{P}\mathrm{O}}_{4 \left(\mathrm{l}\right)}^{3-}\leftrightarrow {\mathrm{K}\mathrm{M}\mathrm{g}\mathrm{P}\mathrm{O}}_{4 \left(\mathrm{s}\right)}$. All these precipitation/dissolution reactions are used to formulate Equation (5), where the precipitation rates of ionic forms of TMs are determined. All already mentioned precipitation/dissolution reactions represent the basis for the determination of the concentrations of precipitates.

The concentrations of the precipitates, formed with ionic forms of appropriate TMs during the batch AD process, are computed as follows:

$${\displaystyle \frac{d{c}_{{i\mathrm{C}\mathrm{O}}_3}}{dt}}=k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},{i\mathrm{C}\mathrm{O}}_3}{\left({\left(c_{i^{2+}}{c}_{{\mathrm{C}\mathrm{O}}_3^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(k_{\mathrm{s}\mathrm{p},{i\mathrm{C}\mathrm{O}}_3}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{M}\mathrm{g},\mathrm{N}\mathrm{i},\mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}\right\}$$

(6)

$${\displaystyle \frac{d{c}_{i_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2}}{dt}}=k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2}{\left({\left({\left(c_{i^{2+}}\right)}^3{\left(c_{{\mathrm{P}\mathrm{O}}_4^{3-}}\right)}^2\right)}^{{\displaystyle \frac{1}{5}}}-{\left(k_{\mathrm{s}\mathrm{p},i_3{\left({\mathrm{P}\mathrm{O}}_4\right)}_2}\right)}^{{\displaystyle \frac{1}{5}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{F}\mathrm{e},\mathrm{N}\mathrm{i}\right\}$$

(7)

$${\displaystyle \frac{d{c}_{i\mathrm{S}}}{dt}}=k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},i\mathrm{S}}{\left({\left(c_{i^{2+}}{c}_{{\mathrm{S}}^{2-}}\right)}^{{\displaystyle \frac{1}{2}}}-{\left(k_{\mathrm{s}\mathrm{p},i\mathrm{S}}\right)}^{{\displaystyle \frac{1}{2}}}\right)}^2, i\in \left\{\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{N}\mathrm{i},\mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}\right\}$$

(8)

The upgraded BioModel, which comprises 80 ordinary differential equations (ODEs) and 54 algebraic equations (AEs), contains 187 model parameters, which are unknown or hard to determine. These model parameters, $x_{\mathrm{m},i}, i=1\dots N_{\mathrm{m}} \left(N_{\mathrm{m}}=187\right)$, which have to be calibrated, can be summarized as follows:

• Initial influent concentrations of 13 types of bacteria in the feeding substrate, $X_j \left(\mathrm{g} {\mathrm{L}}^{-1}\right), j\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$ and concentrations of three types of enzymes in the feeding substrate for degrading carbohydrates, proteins, and lipids, $c_{\mathrm{E}0,i}\left(\mathrm{g} {\mathrm{L}}^{-1}\right), i\in I_{\mathrm{e}\mathrm{n}\mathrm{z}}$;
• Thirteen maximal microbial grow rates at optimal temperature, ${\mu }_{j,\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}}\left({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right), j{\in I}_{\mathrm{b}\mathrm{a}\mathrm{c}};$ 13 microbial decays as a percentage of maximal microbial growth rates $b_{j,\mathrm{d}\mathrm{e}\mathrm{c}} \left(/\right), j{\in I}_{\mathrm{b}\mathrm{a}\mathrm{c}}$; 26 parameters (${\mathrm{p}\mathrm{K}}_j^{\mathrm{l}\mathrm{o}} \left(/\right)$, ${\mathrm{p}\mathrm{K}}_j^{\mathrm{u}\mathrm{p}} \left(/\right)$,$j\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$) which are included in the $\mathrm{p}\mathrm{H}$ function $f_{\mathrm{p}\mathrm{H},j}$ describing the $\mathrm{p}\mathrm{H}$ effects on the growth rate using the Michaelis $\mathrm{p}\mathrm{H}$ function; 13 parameters ${\alpha }_j\left({\mathrm{℃}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right), j\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$; 13 optimal temperatures $T_{j,\mathrm{o}\mathrm{p}\mathrm{t}}\left(℃\right);$ and 13 maximal temperatures $T_{j,\mathrm{m}\mathrm{a}\mathrm{x}}\left(℃\right), j\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$ for the growth rate of $j^{\mathrm{t}\mathrm{h}}$ microbial type;
• Hydrolysis rate constants $k_{\mathrm{h}\mathrm{y}\mathrm{d},i}\left({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right), i{\in I}_{\mathrm{h}\mathrm{y}\mathrm{d}}$ and three Michaelis–Menten half-saturation constants $K_{\mathrm{M}\mathrm{E}0,i} \left(\mathrm{g} {\mathrm{L}}^{-1}\right), i{\in I}_{\mathrm{h}\mathrm{y}\mathrm{d}}$;
• Sixteen Monod saturation constants $k_{\mathrm{M}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ related to various substrates and bacteria: $k_{\mathrm{M},\mathrm{s}\mathrm{u}\mathrm{A}\mathrm{s}\mathrm{u}},$ $k_{\mathrm{M},\mathrm{a}\mathrm{a}\mathrm{A}\mathrm{a}\mathrm{a}},$ $k_{\mathrm{M},\mathrm{g}\mathrm{l}\mathrm{y}\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}},$ $k_{\mathrm{M},\mathrm{o}\mathrm{a}\mathrm{A}\mathrm{o}\mathrm{a}},$ $k_{\mathrm{M},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}},$ $k_{\mathrm{M},\mathrm{b}\mathrm{u}\mathrm{A}\mathrm{b}\mathrm{u}},$ $k_{\mathrm{M},\mathrm{v}\mathrm{a}\mathrm{A}\mathrm{v}\mathrm{a}},$ $k_{\mathrm{M},{\mathrm{H}}_2\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}},$ $k_{\mathrm{M},\mathrm{a}\mathrm{c}\mathrm{M}\mathrm{a}\mathrm{c}},$ $k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{s}},$ $k_{\mathrm{M},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}},$ $k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}},$ $k_{\mathrm{M},\mathrm{a}\mathrm{c}\mathrm{S}\mathrm{a}\mathrm{c}},$ $k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{a}\mathrm{c}},$ $k_{\mathrm{M},{\mathrm{H}}_2\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}},$ and $k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}}$;
• Twenty-six inhibition constants $K_{\mathrm{i}\mathrm{n}\mathrm{h}}\left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ related to (a) VFA inhibition of hydrolysis process $K_{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{V}\mathrm{F}\mathrm{A}}$ and (b) compound and metal ion inhibition of the growth of various bacteria $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2,\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2,\mathrm{A}\mathrm{o}\mathrm{a}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2,\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2,\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2,\mathrm{A}\mathrm{v}\mathrm{a}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{v}\mathrm{a}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{N}\mathrm{H}}_3,\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{C}\mathrm{u}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{Z}\mathrm{n}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{C}\mathrm{r}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{P}\mathrm{b}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{N}\mathrm{i}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{C}\mathrm{u}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{Z}\mathrm{n}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{C}\mathrm{r}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{P}\mathrm{b}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{N}\mathrm{i}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}},$ $K_{\mathrm{i}\mathrm{n}\mathrm{h}{,\mathrm{H}}_2\mathrm{S},\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}},$ $K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{s}},$ $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}},$ $K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{a}\mathrm{c}},$ and $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}},$ and two limitation factors $K_{\mathrm{M}}\left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ of inorganic nitrogen ($K_{\mathrm{M},{\mathrm{N}}_{\mathrm{i}\mathrm{o}}}$) and inorganic phosphorus ($K_{\mathrm{M},{\mathrm{P}}_{\mathrm{i}\mathrm{o}}}$) related to all microbial growth rates;
• Ten parameters of mass transfer rates from liquid to gas phase ${\left(K_{\mathrm{L}}a\right)}_{j,\mathrm{a}}\left({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$ and ${\left(K_{\mathrm{L}}a\right)}_{j,\mathrm{b}}\left({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right), j\in I_{\mathrm{g}\mathrm{a}\mathrm{s}}$;
• Twenty precipitation rate constants $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t}} \left({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right):$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{a}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{u}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{M}\mathrm{g}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{P}\mathrm{b}\mathrm{C}\mathrm{O}}_3},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{Z}\mathrm{n}\mathrm{C}\mathrm{O}}_3},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}\mathrm{S}},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{u}\mathrm{S}},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}\mathrm{S}},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}\mathrm{S}},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{P}\mathrm{b}\mathrm{S}},$ $k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{Z}\mathrm{n}\mathrm{S}},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{a}}_3{\left(\mathrm{P}{\mathrm{O}}_4\right)}_2},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}}_3{\left(\mathrm{P}{\mathrm{O}}_4\right)}_2},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}}_3{\left(\mathrm{P}{\mathrm{O}}_4\right)}_2},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}}_3{\left(\mathrm{P}{\mathrm{O}}_4\right)}_2},$ $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{M}\mathrm{g}\mathrm{N}{\mathrm{H}}_4\mathrm{P}\mathrm{O}}_4}$, and $k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{K}\mathrm{M}\mathrm{g}\mathrm{P}\mathrm{O}}_4}$.

The values of the 20 precipitation rate constants and 10 inhibition constants, which are directly connected to the modeling of TM activities during the AD process, as well as the other 157 model parameters, are determined using a self-developed calibration procedure.

Sensitivity analysis. The sensitivity analysis was performed to determine the so-called importance factor $f_{\mathrm{I}\mathrm{M},i}$ of each $i^{\mathrm{t}\mathrm{h}}$ model parameter from the set of all $N_{\mathrm{m}}$ parameters. For this reason, random values of $5·N_{\mathrm{m}}$ AD model designs were generated (each design is a complete set of all $N_{\mathrm{m}}$ model parameters). The importance of model parameters is determined by their effects on the ${\mathrm{C}\mathrm{H}}_4$ flow rate.

Calibration procedure. In order to calibrate the $N_{\mathrm{m}}$ model parameters, the ASO procedure was used [28,29,30]. The ASO procedure incorporates the upgraded BioModel and the sensitivity analysis with a gradient-based optimization algorithm. The initial values of all model parameters (which are linearly mapped to design variables) were taken as recommended from the available literature. At the beginning of the first stage, some initial and a relatively high activation threshold value $f_T$ was chosen to determine a relatively low number of active design parameters. Within each stage, the active design variables $x_{\mathrm{m},i}^{*}$, for which it holds that the importance factor is greater than the threshold $f_{\mathrm{I}\mathrm{M},i}\ge f_T$, are determined, while all other design variables are designated as passive in the current stage. The values of active design parameters are optimized while keeping the passive ones constant at their recommended values. After that, a new stage with a lower value of threshold $f_T$ is started until all design variables are active and calibrated using the optimization procedure.

BioModel parameter calibration. A unified optimization procedure was used to calibrate the model parameters $x_{\mathrm{m}\mathrm{p},i}, i=1\dots N_{\mathrm{m}}$. The design variables $x_{\mathrm{m},i}$ are defined as normalized model parameters by the relation $x_{\mathrm{m},i}={\displaystyle \frac{x_{\mathrm{m}\mathrm{p},i}-x_{\mathrm{m}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}}{x_{\mathrm{m}\mathrm{p},i}^{\mathrm{U}\mathrm{P}}-x_{\mathrm{m}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}}}$, where $x_{\mathrm{m}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}$ and $x_{\mathrm{m}\mathrm{p},i}^{\mathrm{U}\mathrm{P}}$ are the corresponding lower and upper limits of the $i^{\mathrm{t}\mathrm{h}}$ model parameter $x_{\mathrm{m}\mathrm{p},i}$ to be optimized, respectively. The design variables $x_{\mathrm{m},i}$, the objective function $g_0$, as well as the constraint $g_i$ are defined in various ways, as follows.

The objective function $g_{0,\mathrm{c}\mathrm{a}\mathrm{l}}$ is defined as an integral of daily deviations between the simulated and measured ${\mathrm{C}\mathrm{H}}_4$ flow rates, as shown in Equation (9). These deviations should be minimized. Furthermore, the constraints, which are defined in the standard form $g_{i,\mathrm{c}\mathrm{a}\mathrm{l}}\le 0$, are related to the maximal bacteria concentration by limiting the concentration to below the maximal value of $x_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$, as shown in Equation (10), and to various AD performances by limiting the lower [Equation (11)] and upper [Equation (12)] values of ${\mathrm{C}\mathrm{H}}_4$ flow rates, by limiting the minimal value of the produced ${\mathrm{C}\mathrm{H}}_4$ [Equation (13)], and by prescribing the maximal allowed fractions of the produced ${\mathrm{H}}_2,{\mathrm{H}}_2\mathrm{S},\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{N}{\mathrm{H}}_3$ in biogas [Equations (14)–(16)]. These functions are as described as follows:

$$g_{0,\mathrm{c}\mathrm{a}\mathrm{l}}={\psi }_{0,\mathrm{c}\mathrm{a}\mathrm{l}}\underset{0}{\overset{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\int }}{\left({\displaystyle \frac{Q_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)-Q_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(t\right)}{{\overline{Q}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}}\right)}^{2}dt$$

(9)

$$g_{1,\mathrm{c}\mathrm{a}\mathrm{l}}={\displaystyle \frac{\sum _j\left(X_j-X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\right)}{X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}}, j\in I_{\mathrm{b}\mathrm{a}\mathrm{c}}$$

(10)

$$g_{2,\mathrm{c}\mathrm{a}\mathrm{l}}={\displaystyle \frac{{\int }_0^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\left(0.5+{\tan }^{-1}\left(\frac{10\left(k^{\mathrm{L}\mathrm{O}}{Q}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(t\right)-Q_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)\right)}{{\overline{Q}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}\right)\left(k^{\mathrm{L}\mathrm{O}}{Q}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(t\right)-Q_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)\right)\right)dt}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} {\overline{Q}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}}$$

(11)

$$g_{3,\mathrm{c}\mathrm{a}\mathrm{l}}={\displaystyle \frac{{\int }_0^{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}\left(0.5+{\tan }^{-1}\left(\frac{10\left(Q_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)-k^{\mathrm{U}\mathrm{P}}{Q}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(t\right)\right)}{{\overline{Q}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}\right)\left(Q_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)-k^{\mathrm{U}\mathrm{P}}{Q}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(t\right)\right)\right)dt}{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} {\overline{Q}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}}$$

(12)

$$g_{4,\mathrm{c}\mathrm{a}\mathrm{l}}={\psi }_{4,\mathrm{c}\mathrm{a}\mathrm{l}}\left({\displaystyle \frac{\left|V_{{\mathrm{C}\mathrm{H}}_4}-V_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\right|}{V_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}}}-{\phi }_{{\mathrm{C}\mathrm{H}}_4}\right)$$

(13)

$$g_{5,\mathrm{c}\mathrm{a}\mathrm{l}}={\psi }_{5,\mathrm{c}\mathrm{a}\mathrm{l}}\left({\displaystyle \frac{V_{{\mathrm{H}}_2}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{{\mathrm{H}}_2\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(14)

$$g_{6,\mathrm{c}\mathrm{a}\mathrm{l}}={\psi }_{6,\mathrm{c}\mathrm{a}\mathrm{l}}\left({\displaystyle \frac{V_{{\mathrm{H}}_2\mathrm{S}}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(15)

$$g_{7,\mathrm{c}\mathrm{a}\mathrm{l}}={\psi }_{7,\mathrm{c}\mathrm{a}\mathrm{l}}\left({\displaystyle \frac{V_{{\mathrm{N}\mathrm{H}}_3}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{{\mathrm{N}\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(16)

where ${\psi }_{0,\mathrm{c}\mathrm{a}\mathrm{l}}, {\psi }_{4,\mathrm{c}\mathrm{a}\mathrm{l}}$, ${\psi }_{5,\mathrm{c}\mathrm{a}\mathrm{l}}$, ${\psi }_{6,\mathrm{c}\mathrm{a}\mathrm{l}}$, and ${\psi }_{7,\mathrm{c}\mathrm{a}\mathrm{l}}$ are normalization constants; $t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ is the AD process duration; and the symbols ${\mathrm{Q}}_{{\mathrm{C}\mathrm{H}}_4}\left(t\right) \left(\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$, ${\mathrm{Q}}_{{\mathrm{C}\mathrm{H}}_4,\mathrm{e}\mathrm{x}\mathrm{p}}\left(t\right) \left(\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$, and ${\overline{\mathrm{Q}}}_{{\mathrm{C}\mathrm{H}}_{4,\mathrm{e}\mathrm{x}\mathrm{p}}}\left(\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1}\right)$ denote time-dependent simulated values, time-dependent measured values, and the average values of the measured ${\mathrm{C}\mathrm{H}}_4$ flow rates, respectively. The symbol $X_j \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ denotes the concentration of $j^{\mathrm{t}\mathrm{h}}$ bacteria type from the set $I_{\mathrm{b}\mathrm{a}\mathrm{c}}$, while $X_{\mathrm{b}\mathrm{a}\mathrm{c}}^{\mathrm{m}\mathrm{a}\mathrm{x}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$ denotes the maximal allowed initial concentration of all bacteria. The symbols ${\mathrm{V}}_{{\mathrm{C}\mathrm{H}}_4} \left(\mathrm{L}\right)$ and ${\mathrm{V}}_{{\mathrm{C}\mathrm{H}}_4,\mathrm{e}\mathrm{x}\mathrm{p}} \left(\mathrm{L}\right)$ are the simulated and measured volumes of the produced ${\mathrm{C}\mathrm{H}}_4$, while ${\mathrm{V}}_{{\mathrm{H}}_2}$, ${\mathrm{V}}_{{\mathrm{H}}_2\mathrm{S}},$ ${\mathrm{V}}_{{\mathrm{N}\mathrm{H}}_3}$, and $V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}$ represent the simulated volumes of hydrogen, hydrogen sulfide, ammonia, and biogas, respectively. Finally, the symbol ${\phi }_{{\mathrm{C}\mathrm{H}}_4}$ denotes the prescribed fraction of ${\mathrm{C}\mathrm{H}}_4$, while the symbols ${\phi }_{{\mathrm{H}}_2\mathrm{m}\mathrm{a}\mathrm{x}}$, ${\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}$, and ${\phi }_{{\mathrm{N}\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}$ are the prescribed maximal allowed fractions of hydrogen, hydrogen sulfide, and ammonia, respectively, in the produced biogas. The inverse tangent function in Equations (11) and (12) was employed as a smooth substitute for the conventional step function. This choice is essential to maintain the differentiability of the involved functions, a requirement for using a gradient-based optimizer.

Statistical evaluation methodology. For the evaluation of the proposed upgraded BioModel, the measured AD performance is compared with NS using four statistical indicators (SIs): (i) mean absolute error ${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}}$, (ii) root mean square error ${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$, (iii) coefficient of determination $R^2$, and (iv) the relative index of agreement $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$ [29].

2.3. Optimization of the AD Process

The design variables $x_{\mathrm{p},i}$ were related to some process parameters $x_{\mathrm{p}\mathrm{p},i}, i=1\dots N_{\mathrm{p}}$, and more precisely to the concentrations of TMs ($\mathrm{C}\mathrm{a},\mathrm{C}\mathrm{o},\mathrm{C}\mathrm{r},\mathrm{C}\mathrm{u},\mathrm{F}\mathrm{e},\mathrm{K},\mathrm{M}\mathrm{g},\mathrm{N}\mathrm{a},\mathrm{N}\mathrm{i}, \mathrm{P}\mathrm{b},\mathrm{Z}\mathrm{n}$) in the feeding substrate. Thus, there are a total of 11 design variables $x_{\mathrm{p},i}, i=1\dots N_{\mathrm{p}}, \left(N_{\mathrm{p}}=11\right)$ that are defined as normalized process parameters based on the relation $x_{\mathrm{p},i}={\displaystyle \frac{x_{\mathrm{p}\mathrm{p},i}-x_{\mathrm{p}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}}{x_{\mathrm{p}\mathrm{p},i}^{\mathrm{U}\mathrm{P}}-x_{\mathrm{p}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}}}$, where $x_{\mathrm{p}\mathrm{p},i}^{\mathrm{L}\mathrm{O}}$ and $x_{\mathrm{p}\mathrm{p},i}^{\mathrm{U}\mathrm{P}}$ are the corresponding lower and upper limits of the $i^{\mathrm{t}\mathrm{h}}$ process parameter $x_{\mathrm{p}\mathrm{p},i}, i=1\dots N_{\mathrm{p}}$ to be optimized, respectively.

Three different optimization tasks were formulated using various definitions of objective $g_0$ and constraint $g_i$ functions; these optimization tasks are labeled as Cases A, B, and C. In Case A, a single objective function is defined, while in Cases B and C the multi-objective optimization problems are defined by combining multiple objectives into a single-objective function. This is achieved by summing individual objective functions, multiplied by adequate weighting factors.

In Case A, the produced biogas volume is maximized in the shortest production time possible, as noted in Equation (17). This is achieved using the time-integrated biogas volume, where maximizing the objective function will increase both the total biogas volume and the biogas production rate. In Case B, the produced biogas volume is maximized in the shortest production time possible, and the ${\mathrm{H}}_2\mathrm{S}$ production is minimized, as shown in Equation (18). The objective function, defined by Equation (19), promotes maximal ${\mathrm{C}\mathrm{H}}_4$ volume in the shortest production time possible and minimal ${\mathrm{H}}_2\mathrm{S}$ production. The imposed constraints in Case A are related only to the upper and lower values of the design variables $x_{\mathrm{p},i}, i=1\dots N_{\mathrm{p}}$. In addition to these, the constraint in Case B is related to the maximal allowed ${\mathrm{H}}_2\mathrm{S}$ fraction in biogas, as shown in Equation (20); in Case C, the constraints are related to the maximal allowed fractions of ${\mathrm{H}}_2\mathrm{S}$, ${\mathrm{H}}_2$, and $\mathrm{N}{\mathrm{H}}_3$ in the produced biogas, as shown in Equations (20)–(22). These objective and constraint functions are provided below:

$$g_0=-{\psi }_0\underset{0}{\overset{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\int }}{V}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)dt$$

(17)

$$g_0=-{\psi }_{01}\underset{0}{\overset{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\int }}{V}_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}\left(t\right)+{\psi }_{02}{V}_{{\mathrm{H}}_2\mathrm{S}}$$

(18)

$$g_0=-{\psi }_{01}\underset{0}{\overset{t_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\int }}{V}_{{\mathrm{C}\mathrm{H}}_4}\left(t\right)dt+{\psi }_{02}{V}_{{\mathrm{H}}_2\mathrm{S}}$$

(19)

$$g_1={\psi }_1\left({\displaystyle \frac{V_{{\mathrm{H}}_2\mathrm{S}}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(20)

$$g_2={\psi }_2\left({\displaystyle \frac{V_{{\mathrm{H}}_2}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{{\mathrm{H}}_2,\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(21)

$$g_3={\psi }_3\left({\displaystyle \frac{V_{\mathrm{N}{\mathrm{H}}_3}}{V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}}}}-{\phi }_{\mathrm{N}{\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}\right)$$

(22)

where ${\psi }_0$, ${\psi }_{01}$, ${\psi }_{02}$, ${\psi }_1$, ${\psi }_2$, ${\psi }_3$, and ${\psi }_4$ are the weighting factors; ${\phi }_{{\mathrm{H}}_2\mathrm{S},\mathrm{m}\mathrm{a}\mathrm{x}}$, ${\phi }_{{\mathrm{H}}_2,\mathrm{m}\mathrm{a}\mathrm{x}}$, and ${\phi }_{\mathrm{N}{\mathrm{H}}_3,\mathrm{m}\mathrm{a}\mathrm{x}}$, denote maximal fractions of ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$ in the produced biogas, respectively. The symbols $V_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{a}\mathrm{s}} \left(\mathrm{L}\right)$, $V_{{\mathrm{C}\mathrm{H}}_4} \left(\mathrm{L}\right)$, $V_{{\mathrm{H}}_2} \left(\mathrm{L}\right)$, $V_{{\mathrm{H}}_2\mathrm{S}} \left(\mathrm{L}\right)$, and $V_{\mathrm{N}{\mathrm{H}}_3} \left(\mathrm{L}\right)$ denote the volumes of biogas, ${\mathrm{C}\mathrm{H}}_4$, ${\mathrm{H}}_2$, ${\mathrm{H}}_2\mathrm{S}$, and ${\mathrm{N}\mathrm{H}}_3$, respectively.

2.4. Solution Procedure

The proposed BioModel, the ASO procedure, and the whole optimization procedure were coded in-house in C# language. To solve the system of ODEs, the Runge–Kutta and Euler methods are used. The engaged gradient-based optimization algorithm is based on an approximation method [39,40] that sequentially generates approximate strictly convex and separable nonlinear programming problems and solves them to generate a sequence of converging approximate solutions. The algorithm uses the history of design derivatives of the objective and constraint functions to gradually improve the quality of the approximation. Consequently, the convergence is often relatively fast and stable. Since the analytical derivatives cannot be obtained easily and any derived formulas would be valid only for a particular AD model form, the numerical differentiation using simple forward differences was used in this work to get the needed design derivatives.

All numerical simulations were performed using a desktop computer with Inter i7 3.2 GHz CPU with 8 cores. The CPU time for one simulation of the AD process was approximately 1 s. The CPU time for one full optimization cycle of BioModel parameter calibration, where design derivatives computation was parallelized, is up to 1 min, while 4 s are needed for one optimization cycle of the AD process.

3. Results and Discussion

At first, the results of sensitivity analysis, which were used in BioModel calibration, are presented. Then the AD performance results, obtained by NS using the upgraded BioModel, are compared with experimentally obtained data in the cases of BioModel calibration (bioreactor B1) and BioModel validation (bioreactors B2 and B3) and with various perturbations of calibrated model parameters. Finally, the results of AD model optimization are presented and discussed.

3.1. Sensitivity Analysis

The sensitivity analysis was performed for the AD process in bioreactor B1. The obtained importance factors for 187 BioModel parameters, obtained from sensitivity analysis, are presented in Figure 1.

According to Figure 1, the first ten most important BioModel parameters are the maximal growth rate of methanogenic acetate degraders at optimal temperature ${\mu }_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}}$, lower $\mathrm{p}\mathrm{H}$ drop-off value for methanogenic acetate degraders ${\mathrm{p}\mathrm{K}}_{\mathrm{M}\mathrm{a}\mathrm{c}}^{\mathrm{L}\mathrm{O}}$, inhibition constant for ${\mathrm{H}}_2\mathrm{S}$ by methanogenic acetate degraders $K_{{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{H}}_2\mathrm{S},\mathrm{M}\mathrm{a}\mathrm{c}},$ initial concentration of methanogenic acetate degraders $X_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{i}\mathrm{n}}$, maximal growth rate of acidogenic amino acid degraders at optimal temperature ${\mu }_{\mathrm{A}\mathrm{a}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}}$, inhibition constant for $\mathrm{Z}\mathrm{n}$ by methanogenic acetate degraders $K_{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{Z}\mathrm{n},\mathrm{M}\mathrm{a}\mathrm{c}}$, Monod constant related to acetate and methanogenic acetate degraders $k_{\mathrm{M},\mathrm{a}\mathrm{c}\mathrm{M}\mathrm{a}\mathrm{c}}$, Monod constant related to sulfates and sulfate-reducing bacteria $k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{a}\mathrm{c}}$, maximal growth rate of sulfate-reducing bacteria of acetate at optimal temperature ${\mu }_{\mathrm{S}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}}$, and maximal growth rate of acidogenic sugar degraders at optimal temperature ${\mu }_{\mathrm{A}\mathrm{s}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}}$.

3.2. BioModel Calibration

The calibration of BioModel parameters was performed for the AD process in bioreactor B1. At the beginning, the mean value of interval $\left[x_{\mathrm{m}\mathrm{p},i}^{\mathrm{L}\mathrm{O}},x_{\mathrm{m}\mathrm{p},i}^{\mathrm{U}\mathrm{P}}\right]$ is set as the initial value (initial design) of $i^{\mathrm{t}\mathrm{h}}$ model parameter. The lower, upper, and calibrated values of all 187 BioModel parameters are shown in

Table A1. Lower, upper, and optimal (calibrated) values of BioModel parameters.

$\mathit{i}$	Parameter	${\mathbf{x}}_{\mathbf{m}\mathbf{p},\mathit{i}}^{\mathbf{L}\mathbf{O}}$	${\mathbf{x}}_{\mathbf{m}\mathbf{p},\mathit{i}}^{\mathbf{U}\mathbf{P}}$	Optimal Value	$\mathit{i}$	Parameter	${\mathbf{x}}_{\mathbf{m}\mathbf{p},\mathit{i}}^{\mathbf{L}\mathbf{O}}$	${\mathbf{x}}_{\mathbf{m}\mathbf{p},\mathit{i}}^{\mathbf{U}\mathbf{P}}$	Optimal Value
1	$k_{\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{c}\mathrm{h}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	5.28848	95	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{P}\mathrm{b}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	10	100	55.0005
2	$k_{\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{p}\mathrm{r}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	6.07000	96	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{Z}\mathrm{n}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	10	100	54.9998
3	$k_{\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{l}\mathrm{i}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	6.53405	97	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	500	1000	750.386
4	$K_{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{V}\mathrm{F}\mathrm{A}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.1	0.6	0.27516	98	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}}_3{\left(\mathrm{P}{\mathrm{O}}_{4}4\right)}_2} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	100	74.9785
5	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2,\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}} (\mathrm{g}{\mathrm{L}}^{-1})$	0.01	0.9	0.47075	99	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{a}}_3{\left(\mathrm{P}{\mathrm{O}}_{4}4\right)}_2} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	100	59.3824
6	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2,\mathrm{A}\mathrm{o}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.52991	100	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}}_3{\left(\mathrm{P}{\mathrm{O}}_{4}4\right)}_2} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	100	74.3711
7	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2,\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.56099	101	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}}_3{\left(\mathrm{P}{\mathrm{O}}_{4}4\right)}_2} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	100	74.9990
8	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2,\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.41827	102	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{M}\mathrm{g}\mathrm{N}{\mathrm{H}}_4\mathrm{P}\mathrm{O}}_4} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	200	76.4582
9	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2,\mathrm{A}\mathrm{v}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.44331	103	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{K}\mathrm{M}\mathrm{g}\mathrm{P}\mathrm{O}}_4} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	50	200	124.992
10	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.51747	104	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{s}\mathrm{u}}^{\mathrm{l}\mathrm{o}} (/)$	4.5	5.5	4.98825
11	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.12073	105	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{s}\mathrm{u}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.27958
12	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{A}\mathrm{v}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.09301	106	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{a}\mathrm{a}}^{\mathrm{l}\mathrm{o}} (/)$	4.5	5.5	4.98971
13	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.07857	107	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{a}\mathrm{a}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.48995
14	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.07034	108	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}}^{\mathrm{l}\mathrm{o}} (/)$	4.5	5.5	5.02014
15	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},Ss} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.47785	109	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.22301
16	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.25675	110	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{o}\mathrm{a}}^{\mathrm{l}\mathrm{o}} (/)$	4.5	5.5	5.04878
17	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.31965	111	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{o}\mathrm{a}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.45347
18	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{H}}_2\mathrm{S},\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.53034	112	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.40484
19	$K_{\mathrm{i}\mathrm{n}\mathrm{h},\mathrm{N}{\mathrm{H}}_3,\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.3	0.58207	113	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}}^{\mathrm{u}\mathrm{p}} (/)$	8.0	9.0	8.20391
20	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{C}\mathrm{u}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.41420	114	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{b}\mathrm{u}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.31587
21	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{Z}\mathrm{n}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.35448	115	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{b}\mathrm{u}}^{\mathrm{u}\mathrm{p}} (/)$	8.0	9.0	8.61508
22	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{C}\mathrm{r}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.44256	116	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{v}\mathrm{a}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.13055
23	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{P}\mathrm{b}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.45450	117	${\mathrm{p}\mathrm{K}}_{\mathrm{A}\mathrm{v}\mathrm{a}}^{\mathrm{u}\mathrm{p}} (/)$	8.0	9.0	8.53030
24	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{N}\mathrm{i}}^{2+},\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.45174	118	${\mathrm{p}\mathrm{K}}_{\mathrm{M}\mathrm{a}\mathrm{c}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.00441
25	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{C}\mathrm{u}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.51131	119	${\mathrm{p}\mathrm{K}}_{\mathrm{M}\mathrm{a}\mathrm{c}}^{\mathrm{u}\mathrm{p}} (/)$	8.0	9.0	8.34907
26	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{Z}\mathrm{n}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.58076	120	${\mathrm{p}\mathrm{K}}_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.04491
27	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{C}\mathrm{r}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.43723	121	${\mathrm{p}\mathrm{K}}_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}}^{\mathrm{u}\mathrm{p}} (/)$	8.0	9.0	8.48158
28	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{P}\mathrm{b}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.45489	122	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{s}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	5.96875
29	$K_{\mathrm{i}\mathrm{n}\mathrm{h},{\mathrm{N}\mathrm{i}}^{2+},\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.49740	123	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{s}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.07330
30	$K_{\mathrm{M},{\mathrm{N}}_{\mathrm{i}\mathrm{o}}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.001	0.01	0.00513	124	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.21593
31	$K_{\mathrm{M},{\mathrm{P}}_{\mathrm{i}\mathrm{o}}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.001	0.01	0.00536	125	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	8.00344
32	$k_{\mathrm{M},\mathrm{s}\mathrm{u}\mathrm{A}\mathrm{s}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.45069	126	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{a}\mathrm{c}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	5.64308
33	$k_{\mathrm{M},\mathrm{a}\mathrm{a}\mathrm{A}\mathrm{a}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.38091	127	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{a}\mathrm{c}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	7.64437
34	$k_{\mathrm{M},\mathrm{g}\mathrm{l}\mathrm{y}\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.28036	128	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}}^{\mathrm{l}\mathrm{o}} (/)$	5.5	6.5	6.02774
35	$k_{\mathrm{M},\mathrm{o}\mathrm{a}\mathrm{A}\mathrm{o}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.02559	129	${\mathrm{p}\mathrm{K}}_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}}^{\mathrm{u}\mathrm{p}} (/)$	7.5	8.5	7.98378
36	$k_{\mathrm{M},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.72705	130	${\alpha }_{\mathrm{A}\mathrm{s}\mathrm{u}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
37	$k_{\mathrm{M},\mathrm{b}\mathrm{u}\mathrm{A}\mathrm{b}\mathrm{u}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.86578	131	${\alpha }_{\mathrm{A}\mathrm{a}\mathrm{a}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
38	$k_{\mathrm{M},\mathrm{v}\mathrm{a}\mathrm{A}\mathrm{v}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.65487	132	${\alpha }_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
39	$k_{\mathrm{M},{\mathrm{H}}_2\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.89994	133	${\alpha }_{\mathrm{A}\mathrm{o}\mathrm{a}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
40	$k_{\mathrm{M},\mathrm{a}\mathrm{c}\mathrm{M}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.70222	134	${\alpha }_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
41	$k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{s}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.09066	135	${\alpha }_{\mathrm{A}\mathrm{b}\mathrm{u}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
42	$k_{\mathrm{M},\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.58511	136	${\alpha }_{\mathrm{A}\mathrm{v}\mathrm{a}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
43	$k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.61723	137	${\alpha }_{\mathrm{M}\mathrm{a}\mathrm{c}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00054
44	$k_{\mathrm{M},\mathrm{a}\mathrm{c}\mathrm{S}\mathrm{a}\mathrm{c}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.24452	138	${\alpha }_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00015	0.00019	0.00055
45	$k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{a}\mathrm{c}}$(g L−1)	0.01	0.9	0.24192	139	${\alpha }_{\mathrm{S}\mathrm{s}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
46	${\mathrm{k}}_{\mathrm{M},{\mathrm{H}}_2\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.33351	140	${\alpha }_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
47	$k_{\mathrm{M},\mathrm{S}\mathrm{i}\mathrm{o},\mathrm{a}\mathrm{t}\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01	0.9	0.34444	141	${\alpha }_{\mathrm{S}\mathrm{a}\mathrm{c}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
48	${\mu }_{\mathrm{A}\mathrm{s}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	7.27177	142	${\alpha }_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} ({\mathrm{K}}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.00016	0.00020	0.00055
49	${\mu }_{\mathrm{A}\mathrm{a}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	9.99998	143	$T_{\mathrm{A}\mathrm{s}\mathrm{u},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9854
50	${\mu }_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	9.18755	144	$T_{\mathrm{A}\mathrm{s}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
51	${\mu }_{\mathrm{A}\mathrm{o}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	9.91931	145	$T_{\mathrm{A}\mathrm{s}\mathrm{a},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9685
52	${\mu }_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	3.51354	146	$T_{\mathrm{A}\mathrm{a}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
53	${\mu }_{\mathrm{A}\mathrm{b}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	1.62130	147	$T_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9975
54	${\mu }_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	2.94090	148	$T_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
55	${\mu }_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	5.46362	149	$T_{\mathrm{A}\mathrm{o}\mathrm{a},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9941
56	${\mu }_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	3.64307	150	$T_{\mathrm{A}\mathrm{o}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
57	${\mu }_{\mathrm{S}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	5.63776	151	$T_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9866
58	${\mu }_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	1.17762	152	$T_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
59	${\mu }_{\mathrm{S}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	9.34170	153	$T_{\mathrm{A}\mathrm{b}\mathrm{u},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	55	65	60.0039
60	${\mu }_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x},T_{\mathrm{o}\mathrm{p}\mathrm{t}}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	6.70977	154	$T_{\mathrm{A}\mathrm{b}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	65	75	70.0000
61	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{s}\mathrm{u}} (/)$	0.01	0.05	0.02996	155	$T_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	55	65	60.0130
62	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{a}\mathrm{a}} (/)$	0.01	0.05	0.02684	156	$T_{\mathrm{A}\mathrm{v}\mathrm{a},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	65	75	70.0000
63	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}} (/)$	0.01	0.05	0.02985	157	$T_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	54.9936
64	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{o}\mathrm{a}} (/)$	0.01	0.05	0.03009	158	$T_{\mathrm{M}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
65	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} (/)$	0.01	0.05	0.01590	159	$T_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	50	60	55.0023
66	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{b}\mathrm{u}} (/)$	0.01	0.05	0.04289	160	$T_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0000
67	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{A}\mathrm{v}\mathrm{a}} (/)$	0.01	0.05	0.02859	161	$T_{\mathrm{S}\mathrm{s},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	30	40	35.1549
68	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{M}\mathrm{a}\mathrm{c}} (/)$	0.01	0.05	0.01980	162	$T_{\mathrm{S}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0022
69	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}} (/)$	0.01	0.05	0.04300	163	$T_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	30	40	35.0021
70	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{S}\mathrm{s}} (/)$	0.01	0.05	0.03071	164	$T_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o},\mathrm{m}\mathrm{a}\mathrm{x}}$(°C)	60	70	65.0000
71	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} (/)$	0.01	0.05	0.03000	165	$T_{\mathrm{S}\mathrm{a}\mathrm{c},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	30	40	36.2293
72	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{S}\mathrm{a}\mathrm{c}} (/)$	0.01	0.05	0.02867	166	$T_{\mathrm{S}\mathrm{a}\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	65.0008
73	$b_{\mathrm{d}\mathrm{e}\mathrm{c},\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} (/)$	0.01	0.05	0.02999	167	$T_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{o}\mathrm{p}\mathrm{t}} (℃)$	30	40	34.9999
74	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{C}\mathrm{O}}_2,\mathrm{a}} ({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	5	3.59248	168	$T_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x}} (℃)$	60	70	64.9999
75	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{C}\mathrm{O}}_2,\mathrm{b}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	10	20	14.97329	169	$K_{\mathrm{M}\mathrm{E}0,\mathrm{c}\mathrm{h}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00009
76	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{C}\mathrm{H}}_4,\mathrm{a}} ({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	5	2.97810	170	$K_{\mathrm{M}\mathrm{E}0,\mathrm{p}\mathrm{r}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00024
77	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{C}\mathrm{H}}_4,\mathrm{b}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	10	20	14.99536	171	$K_{\mathrm{M}\mathrm{E}0,\mathrm{l}\mathrm{i}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00053
78	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{H}}_2,\mathrm{a}} ({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0001	0.001	0.00086	172	$X_{\mathrm{A}\mathrm{s}\mathrm{u}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.26646
79	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{H}}_2,\mathrm{b}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.001	0.01	0.00641	173	$X_{\mathrm{A}\mathrm{a}\mathrm{a}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.49999
80	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{H}}_2\mathrm{S},\mathrm{a}} ({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0001	0.001	0.00053	174	$X_{\mathrm{A}\mathrm{g}\mathrm{l}\mathrm{y}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.35606
81	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{H}}_2\mathrm{S},\mathrm{b}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.001	0.01	0.00565	175	$X_{\mathrm{A}\mathrm{o}\mathrm{a}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.48358
82	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{N}\mathrm{H}}_3,\mathrm{a}} ({℃}^{-1} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0001	0.001	0.00055	176	$X_{\mathrm{A}\mathrm{p}\mathrm{r}\mathrm{o}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.10110
83	${\left(K_{\mathrm{L}}a\right)}_{{\mathrm{N}\mathrm{H}}_3,\mathrm{b}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.001	0.01	0.00550	177	$X_{\mathrm{A}\mathrm{b}\mathrm{u}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.20184
84	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{a}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	2.26219	178	$X_{\mathrm{A}\mathrm{v}\mathrm{a}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.14735
85	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{M}\mathrm{g}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	6.23296	179	$X_{\mathrm{M}\mathrm{a}\mathrm{c}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.10001
86	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	2.45088	180	$X_{\mathrm{M}\mathrm{h}\mathrm{y}\mathrm{d}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.31548
87	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	5.34170	181	$X_{\mathrm{S}\mathrm{p}\mathrm{r}\mathrm{o}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.25498
88	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{u}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	2.21574	182	$X_{\mathrm{S}\mathrm{a}\mathrm{c}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.10000
89	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{P}\mathrm{b}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	5.47297	183	$X_{\mathrm{S}\mathrm{s}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.41458
90	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{Z}\mathrm{n}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	2.61502	184	$X_{\mathrm{S}\mathrm{h}\mathrm{y}\mathrm{d}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.1	0.5	0.34090
91	$k_{{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{o}\mathrm{C}\mathrm{O}}_3} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1	10	4.96625	185	$c_{\mathrm{E}0,\mathrm{c}\mathrm{h}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00063
92	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{F}\mathrm{e}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	1000	5000	2999.88	186	$c_{\mathrm{E}0,\mathrm{p}\mathrm{r}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00046
93	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{N}\mathrm{i}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	500	1000	2999.88	187	$c_{\mathrm{E}0,\mathrm{l}\mathrm{i}} \left(\mathrm{g} {\mathrm{L}}^{-1}\right)$	0.00001	0.001	0.00054
94	$k_{\mathrm{c}\mathrm{r}\mathrm{y}\mathrm{s}\mathrm{t},\mathrm{C}\mathrm{u}\mathrm{S}} ({\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	500	1000	749.971

in Appendix A.

For the selection of active design variables, the threshold values of $f_T$ were sequentially chosen as 0.2 (Set 1), 0.1 (Set 2), 0.001 (Set 3), and 0 (Set 4, optimal design). For each set, the total number $N_{{x_{\mathrm{m},i}}^{*}}$ of the corresponding active design variables $x_{\mathrm{m},i}^{*}$, number of iterations $N_{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}}$ during the calibration, and the value of objective function $g_0$ are collected in

Table 3. Active sets data related to the ASO procedure.

Set	${\mathit{f}}_{\mathit{T}}$	${\mathit{N}}_{{{\mathit{x}}_{\mathbf{m},\mathit{i}}}^{\mathit{*}}}$	${\mathit{N}}_{\mathbf{i}\mathbf{t}\mathbf{e}\mathbf{r}}$	${\mathit{g}}_0$
1	0.2	17	16	0.25700
2	0.1	30	59	0.05285
3	0.001	89	78	0.02280
4 (ODcal)	0	187	190	0.01185

. The optimized values of active design variables $x_{\mathrm{m},i}^{*}$, obtained by considering Set 1, were used as initial values of these variables for the optimization of the next set. The value of the objective function was gradually minimized from the initial design (ID_cal) value of 3.1060 to the optimal design (OD_cal, Set 4) value of 0.01185.

The results of BioModel calibration by the ASO procedure are presented in Figure 2. It can be seen that the differences between measured and simulated ${\mathrm{C}\mathrm{H}}_4$ flow rates in the initial design are extremely high. Furthermore, the simulated values are becoming closer to experimental data as the number of active parameters increases. This is confirmed also by the average difference between simulated and measured ${\mathrm{C}\mathrm{H}}_4$ flow rates (average absolute daily difference divided by average measurement value). The average difference of ${\mathrm{C}\mathrm{H}}_4$ flow rates is reduced during calibration from 95.5% (initial values of model parameters; ID_cal), through 40.9% (Set 1), 23.6% (Set 2), and 16.8% (Set 3) to 12.0% (final calibrated values of model parameters; OD_cal). It is clearly evident that the simulated ${\mathrm{C}\mathrm{H}}_4$ flow rates, obtained under the optimal design with all calibrated BioModel parameters, agree very well with the measured values throughout the whole time period of the AD process.

The upgraded BioModel enables predictions of the concentration of TEs during the AD process. As an example, the concentrations of various Fe and Mg species at OD_cal are presented in Figure 3.

The calibrated values of BioModel parameters are confirmed also by error, ${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}}$, and ${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$ as well as efficiency, $R^2$ and $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$, SIs, as shown in

Table 4. Statistical indicators of CH₄ flow rate; bioreactor B1, BioModel calibration.

SIs	BioModel Calibration, Bioreactor B1
	IDcal	Set 1	Set 2	Set 3	ODcal
${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}} (\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0583	0.0250	0.0144	0.0102	0.0074
${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}} (\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.2933	0.08438	0.0382	0.0252	0.0181
$R^2 (/)$	0.2374	0.6491	0.9360	0.9810	0.9855
$I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}} (/)$	0.4872	0.8440	0.9774	0.9902	0.9953

. For bioreactor B1 (BioModel calibration), the calculated values of SIs are related to methane flow rate $Q_{{\mathrm{C}\mathrm{H}}_4}$. During calibration process using the ASO procedure, from the ID_cal to the OD_cal design, the values of ${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}}$ and ${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$ decrease by approximately eightfold and sixteenfold, respectively. On the other hand, the values of $R^2$ and $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$ increase by approximately fourfold and twofold, respectively. The obtained values of all error and efficiency SIs for OD_cal design are acceptable.

3.3. BioModel Validation

The validation of BioModel parameters, which were calibrated with respect to the experimental data of the AD process in bioreactor B1, is performed for the AD processes in bioreactors B2 and B3. For these two cases, the numerically obtained biogas flow rates are compared with measured data, as shown in Figure 4. The average differences in biogas flow rates are 11.5% (bioreactor B2) and 11.2% (bioreactor B3). In both cases, the biogas flow rate obtained using NS agree very well with the experimentally obtained one. This good agreement of AD performance confirms also the efficiency of the ASO procedure for calibration.

Furthermore, with respect to the measured and simulated biogas flow rates in bioreactors B2 and B3 (BioModel validation), various SIs were calculated, as shown in

Table 5. Statistical indicators of biogas flow rate; bioreactors B2 and B3, BioModel validation.

SIs	B2	B3
${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}} (\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0033	0.0654
${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}} (\mathrm{L} {\mathrm{d}\mathrm{a}\mathrm{y}}^{-1})$	0.0111	0.1817
$R^2 (/)$	0.9780	0.9775
$I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}} (/)$	0.9877	0.9931

. All obtained values of the error SIs are acceptable.

3.4. Importance of Accurate Calibration of BioModel Parameters

In order to test the importance of the accuracy of calibration of all model parameters, an extra study was performed to include various perturbations of the calibrated values; the observed bioreactors are B2 and B3. For this purpose, the impact of random perturbations of all BioModel parameters on the AD performance was investigated. In the case of small perturbations, the values of model parameters were randomly generated based on a $\pm 1\%$ perturbation of their calibrated values (OD_cal $\pm 1\%$). In the case of high perturbation, the values of model parameters were randomly generated based on a $\pm 5\%$ perturbation of their calibrated values (OD_cal $\pm 5\%$). The comparison of the simulated AD performances with the measured data are shown in Figure 5. The values of average differences in biogas flow rates in the case of validation in bioreactor B2 increases by increasing the perturbation of calibrated values of model parameters from 11.5% (OD_cal) and 12.0% (OD_cal $\pm 1\%$) to 29.1% (OD_cal $\pm 5\%$). In the case of validation in bioreactor B3, the obtained average difference of the biogas flow rate increases from 11.2% (OD_cal) and 14.6% (OD_cal $\pm 1\%$) to 34.5% (OD_cal $\pm 5\%$). It is evident that using 5% random perturbation of calibrated values of model parameters, the discrepancy between simulated and experimental AD performance can increase up to threefold with respect to those obtained with OD_cal.

The comparison of biogas flow rates in Figure 5 clearly shows that the disagreement increases by increasing the perturbation of model parameters. A similar conclusion can be made also by comparing the SIs, as shown in

Table 6. Statistical indicators of biogas flow rate; perturbation of model parameters, bioreactors B2 and B3, validation period.

Bioreactor	Design	SIs
		${\mathit{\epsilon }}_{\mathbf{M}\mathbf{A}\mathbf{E}} (\mathbf{L} {\mathbf{d}\mathbf{a}\mathbf{y}}^{-1})$	${\mathit{\epsilon }}_{\mathbf{R}\mathbf{M}\mathbf{S}\mathbf{E}} (\mathbf{L} {\mathbf{d}\mathbf{a}\mathbf{y}}^{-1})$	${\mathit{R}}^2 (/)$	${\mathit{I}}_{\mathbf{A},\mathbf{r}\mathbf{e}\mathbf{l}} (/)$
B2	ODcal	0.0033	0.0111	0.9780	0.9877
	ODcal ± 1%	0.0036	0.0113	0.9712	0.9866
	ODcal ± 5%	0.0087	0.0303	0.8681	0.9261
B3	ODcal	0.0654	0.1817	0.9775	0.9931
	ODcal ± 1%	0.0850	0.2544	0.9571	0.9869
	ODcal ± 5%	0.2010	0.6241	0.8090	0.9311

.

Increasing perturbations of the calibrated values of model parameters results in progressively degraded accuracy of simulated performances. The values of ${\epsilon }_{\mathrm{M}\mathrm{A}\mathrm{E}}$ and ${\epsilon }_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$ increase, while the values of $R^2$ and $I_{\mathrm{A},\mathrm{r}\mathrm{e}\mathrm{l}}$ decrease. The results in Figure 5 and Table 6 show that accurate calibration of all model parameters is very important.

3.5. Optimized AD Process

In order to improve AD performance, the concentrations of all considered light TMs (Ca, K, Mg, and Na) and heavy TMs (Co, Cr, Cu, Fe, Ni, Pb, and Zn) are optimized. The prescribed allowed lower and upper values of TMs concentrations, actual initial TMs concentrations in the feeding substrate of bioreactor B1, as well as the optimal values of TMs in optimization Cases A, B, and C are given in

Table 7. TM concentration; lower, upper, initial, and optimized values.

No.	TM	${\mathbf{x}}_{\mathbf{p}\mathbf{p},\mathit{i}}^{\mathbf{L}\mathbf{O}}$	${\mathbf{x}}_{\mathbf{p}\mathbf{p},\mathit{i}}^{\mathbf{U}\mathbf{P}}$	ID	OD Case A	OD Case B	OD Case C
1	$c_{\mathrm{C}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.30000	4.00000	3.00000	3.9870	2.2957	3.1290
2	$c_{\mathrm{K}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.32000	4.20000	3.20000	3.2587	3.2263	3.1430
3	$c_{\mathrm{M}\mathrm{g}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.08000	1.80000	0.80000	1.7868	1.7970	0.0887
4	$c_{\mathrm{N}\mathrm{a}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.03000	0.80000	0.30000	0.3000	0.3000	0.3000
5	$c_{\mathrm{F}\mathrm{e}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.01100	5.11000	0.11000	5.0431	5.1100	5.1099
6	$c_{\mathrm{C}\mathrm{o}} (\mathrm{g} {\mathrm{L}}^{-1})$	0.00200	0.12000	0.02000	0.0251	0.0229	0.0239
7	$c_{\mathrm{C}\mathrm{r}} (\mathrm{m}\mathrm{g} {\mathrm{L}}^{-1})$	0.30000	103.000	3.00000	0.3342	0.3000	0.3000
8	$c_{\mathrm{N}\mathrm{i}} (\mathrm{m}\mathrm{g} {\mathrm{L}}^{-1})$	0.10000	101.000	1.00000	0.1113	0.1000	0.1000
9	$c_{\mathrm{P}\mathrm{b}} (\mathrm{m}\mathrm{g} {\mathrm{L}}^{-1})$	0.02000	100.200	0.20000	0.0218	0.0200	0.0200
10	$c_{\mathrm{C}\mathrm{u}} (\mathrm{m}\mathrm{g} {\mathrm{L}}^{-1})$	1.00000	110.000	10.00000	1.1192	1.0001	1.0001
11	$c_{\mathrm{Z}\mathrm{n}} (\mathrm{m}\mathrm{g} {\mathrm{L}}^{-1})$	5.00000	150.000	50.00000	5.9228	5.0003	5.0010

. The lower values are prescribed as 10% of the initial concentration in the feeding substrate; these TM concentration values can be achieved by removing 90% of TMs from the feeding substrate. The upper values of TM concentrations can be achieved by adding TMs to the feeding substrate. As it is evident from Table 7, only the initial value of Na concentration remains the same in all optimization cases. This result is not surprising because the Na concentration in the feeding substrate in the initial design is already in the optimal range for the growth of mesophilic hydrogenotrophic methanogens and the associated methane production [19].

Some results of AD process optimization are shown in Figure 6, Figure 7, Figure 8 and Figure 9. The optimized amounts of added/removed TMs in the optimization Cases A, B, and C with respect to the initial TM concentrations in the feeding substrate are shown in Figure 6. It is evident that in all cases, the heavy metals Cr, Ni, Pb, Cu, and Zn have to be removed from the feeding substrate before starting the AD process, while Fe and Co should be added. The light TMs should be added or removed independently of the optimization case. All light TMs should be added in all cases, with the exceptions of removing Ca in Case B as well as Mg and K in Case C.

Using these optimized concentrations of TMs in the feeding substrate, the quality and quantity of the produced biogas are improved in all cases, as shown in Figure 7. The obtained results show that the highest biogas production is reached in Case A, where the objective of optimization is related only to the maximization of biogas volume. As expected, by including other requirements in the objective and constraint functions (Cases B and C), the improvement in biogas production becomes smaller. By comparing the obtained optimal designs with the initial design in Figure 8, the biogas production increases by 22.6%, 17.3%, and 5.4% in Cases A, B, and C, respectively. The CH₄ production with respect to the initial design increases by 29.4%, 23.6%, and 13.2% in Cases A, B, and C, respectively.

By analyzing the results presented in Figure 7, it can be concluded that the biogas quality increases from Case A through Case B to Case C. For example, the produced biogas in Case C contains 7.4% more CH₄, 31.9% less H₂, 46.5% less H₂S, and 45.4% less NH₃ with respect to the initial design, as shown in Figure 8. The highest content of CH₄ in Case C is a consequence of lower concentrations of Ca, Mg, and K in the feeding substrate. Namely, the obtained optimal Ca²⁺, Mg²⁺, and K⁺ concentrations are important for dehydrogenase activity, cell viability, and maintenance of the morphology of microbes. Higher K⁺ concentrations could lead to the displacement of exchangeable metal ions from exchangeable sites and consequently to the loss of micronutrients important for methanogens activity, or to cell dehydration, leading to the inhibition of methanogen activity [19]. Since the methanogens are less tolerant to high Mg²⁺ and Ca²⁺ concentrations, the maximization of CH₄ production in Case C results in the lowest optimized Mg²⁺ concentration.

The main impact for the reduction of the content of H₂S in biogas is due to the added amount of Fe²⁺ in the feeding substrate. The Fe²⁺ ions form various precipitates, such as carbonates, sulfates, and phosphates. The concentrations of various iron species during the AD process are presented in Figure 9. It is clearly evident that the formation of iron precipitates (sulfide, carbonate, and phosphate) is the fastest in the initial stage of the AD process in Case C (at optimal design) and that the amounts of FeS formed during the AD process at optimal designs in Cases A, B, and C are significantly larger than those obtained with the initial designs. Consequently, the content of H₂S in the produced biogas is the smallest in Case C; the reduction of H₂S with respect to initial design is 46.5%, followed by 42.5% in Case B and 38.6% in Case A, as shown in Figure 8.

An analysis of the obtained results of the AD process optimization shows the following findings: (i) the amount of H₂S in biogas depends not only on the availability of Fe²⁺ ions but also on the formation of FeS at the initial stage of the AD process; (ii) the necessary reductions in metal concentrations in the feeding substrate, obtained during AD process optimization, require appropriate pretreatments of the feeding substrate, including adsorption of various light and heavy metals by various adsorbents (biochar, active carbon, nanomaterials, etc.) or by the addition of various substrates with lower content of metals; (iii) by maximizing the production of CH₄ and minimizing the content of H₂S simultaneously and by introducing the limits on the content of H₂, H₂S, and NH₃ by the constraint functions (Case C), the optimal concentrations of TMs increased biogas and CH₄ production by 5.4% and 13.5%, respectively, while H₂, H₂S, and NH₃ decreased by 28.2%, 43.6%, and 42.5%, respectively. Furthermore, the content of CH₄ increases by 7.4%, while H₂, H₂S, and NH₃ decreases by 31.9%, 46.5%, and 45.4%, respectively.

4. Conclusions

This study deals with the modeling, calibration, validation, and numerical optimization of the AD process with special attention given to the role of TMs in batch anaerobic digestion. For this reason, a mechanistically inspired BioModel was enhanced by modeling TM activities in biochemical, chemical, and physicochemical processes; all 187 model parameters, including TM-related parameters, were calibrated by considering their interactions. The average differences in the observed AD performance in bioreactors B1, B2, and B3 are $11.6\%\pm 0.4\%$ (calibrated model parameters, OD_cal). A relatively small increase in average differences in biogas flow rate in bioreactors B2 and B3 is obtained at low perturbation of calibrated parameters. Meanwhile, at high perturbation (OD_cal $\pm 5\%$), the average differences increase up to 34.5%; a random ±5% perturbation of calibrated values of model parameters results in rather inaccurate predictions of AD performance. Finally, by optimizing TM concentrations in the feeding substrate, the quantity and quality of the produced biogas can be improved essentially; the biogas and CH₄ production increased by 5.4% and 13.5%, respectively, while H₂, H₂S, and NH₃ decreased by 28.2%, 43.6%, and 42.5%, respectively. In the future, further research focusing on the effects of multiple metals in bioreactors on the AD process as well as on the specific thresholds of the stimulatory and inhibitory concentrations during various stages of the AD process would be highly useful.