# A general framework to model the fate of trace elements in anaerobic digestion environments

Bikash Chandra Maharaj, Maria Rosaria Mattei, Luigi Frunzo, Eric D. van Hullebusch, Giovanni Esposito

► To cite this version:

Bikash Chandra Maharaj, Maria Rosaria Mattei, Luigi Frunzo, Eric D. van Hullebusch, Giovanni Esposito. A general framework to model the fate of trace elements in anaerobic digestion environments. Scientific Reports, 2021, 11, 19 pp. 10.1038/s41598-021-85403-2. insu-03590103

HAL Id: insu-03590103

https://insu.hal.science/insu-03590103

Submitted on 26 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A general framework to model the fate of trace elements in anaerobic digestion environments

Bikash Chandra Maharaj<fnref n="1" />,²,³✉, Maria Rosaria Mattei², Luigi Frunzo², Eric D. van Hullebusch⁴ & Giovanni Esposito³

Due to the multiplicity of biogeochemical processes taking place in anaerobic digestion (AD) systems and limitations of the available analytical techniques, assessing the bioavailability of trace elements (TEs) is challenging. Determination of TE speciation can be facilitated by developing a mathematical model able to consider the physicochemical processes affecting TEs dynamics. A modeling framework based on anaerobic digestion model no 1 (ADM1) has been proposed to predict the biogeochemical fate TEs in AD environments. In particular, the model considers the TE adsorption-desorption reactions with biomass, inerts and mineral precipitates, as well as TE precipitation/dissolution, complexation reactions and biodegradation processes. The developed model was integrated numerically, and numerical simulations have been run to investigate the model behavior. The simulation scenarios predicted the effect of (i) organic matter concentration, (ii) initial TEs concentrations, (iii) initial Ca-Mg concentrations, (iv) initial EDTA concentration, and (v) change in TE binding site density, on cumulative methane production and TE speciation. Finally, experimental data from a real case continuous AD system have been compared to the model predictions. The results prove that this modelling framework can be applied to various AD operations and may also serve as a basis to develop a model-predictive TE dosing strategy.

Anaerobic digestion (AD) of organic wastes is a well-established alternative source of sustainable energy<fnref n="1" />. Efforts have been devoted to improve the biogas production by recruiting process optimization and scale up techniques²⁻⁷. Dosing of trace elements (TEs: Fe, Ni, Co, Se, Mn, Zn and Mo) in AD is one of such avenues⁸⁻¹¹. This is primarily due to the sensitivity of the microbial community to TEs, in particular acetogens and methanogens. Indeed, the syntrophy between acetogens and methanogens is largely attributed to optimum availability and supply of TEs. Any lack or excess of TEs in the reactor leads to acidification and toxicity, respectively. Moreover, the complex biogeochemistry of TEs requires a thorough analysis regarding the fate of TEs in AD systems¹². Examining the biogeochemistry and the dynamic behavior of TEs in AD is necessary to individuate important process parameters which would be the first step towards a model predictive TE dosage control.

TEs are structural components of important enzymes involved in AD biochemical processes¹³. Fe is a structural component of the majority of enzymes in AD pathways¹⁴,¹⁵. Also, it serves as a redox carrier¹⁶. Ni is a part of coenzyme F₄₃₀ ring structure of methyl reductase which facilitates reductions of methyl coenzyme M to methane during methanogenesis¹⁷. Co is essential to both acetoclastic and hydrogenotrophic methanogens. It is an integral part of coronoid structure of vitamin B12 which binds to coenzyme methylase¹⁸. Thus, these bioavailable TEs are necessary micronutrients for the consortium of microorganisms present in AD.

The bioavailable fraction of TEs in an AD system is a function of the bio-uptake, precipitation, complexation and adsorption processes¹⁹,²⁰. Proper bioavailability of TEs in an AD system is necessary for an efficient methane production²¹. Microbial growth depends on the uptake and transport of these essential TEs from the bulk solution¹⁹. Quantification of mineral precipitation/dissolution²²⁻³⁰, organic complexation, and adsorption/desorption is crucial to determine the TE bioavailability in the liquid phase.

<footnote n="1">¹Department of Civil and Mechanical Engineering, University of Cassino and the Southern Lazio, Via di Biasio 43, 03043 Cassino, FR, Italy. ²Department of Mathematics and Applications “Renato Cacciopoli”, University of Naples Federico II, Naples, Italy. ³Department of Civil, Architectural and Environmental Engineering, University of Naples Federico II, Naples, Italy. ⁴Institut de Physique du Globe de Paris, Université de Paris, CNRS, F-75005 Paris, France. ✉email: bikash.maharaj@gmail.com</footnote>
Due to the significant role played by TEs in AD systems, the difficulty in measuring the bioavailable TE fractions<sup>31,32</sup>, as well as the intricate TE biogeochemistry, a mathematical model able to simulate the complex dynamics of TEs in anaerobic environments would be a useful tool<sup>21</sup>. Although there have been experimental efforts to quantify TE in the liquid phase of anaerobic reactors, a mathematical model-based approach will augment and accelerate the TE quantification towards an accurate and reliable dosing strategy designed on case to case basis.

In this regard, incorporating physicochemical processes in AD modelling is crucial. Many studies have highlighted the omission of important physicochemical processes (other than gas transfer and acid base transformation) in ADM1<sup>33</sup>. Subsequently, some attempts have been made to improve the physicochemical predictions of the model. In a minor model improvement calcium carbonate (CaCO<sub>3</sub>) precipitation was included to demonstrate an industrial application of the model<sup>23</sup>. Later on, a more detailed physicochemical framework was implemented to study the precipitation of phosphorous with alkaline earth metals (Ca and Mg)<sup>27</sup>. Similarly, the interactions of phosphorous, sulfur and iron were studied<sup>25</sup> with emphasis on phosphorous recovery and iron precipitation as one of the processes affecting microbial kinetics. Evidently, these physicochemical improvements mainly focused on the reactions of phosphorous metabolism. Moreover, only selected precipitation/dissolution reactions have been considered by these models while neglecting other major physicochemical processes such as complexation and adsorption/desorption. Recently, a model for TE quantification has been proposed which is mechanistic in nature and considers only representative components and species<sup>34</sup>, thus lacking a detailed description of major species and components required to understand TE biogeochemistry in AD. Apart from these, there have been attempts to describe complexation<sup>35</sup> and adsorption processes<sup>36</sup> in wastewater systems but they are not specific to anaerobic environments. Ultimately, these studies would not be sufficient to understand the complex mechanisms associated with TEs in AD as, more often, they lack one or many components and processes necessary to describe the TE biogeochemistry in AD systems<sup>26,37</sup>.

In contrast, numerous scientific reports<sup>20,27,33,38–41</sup> have highlighted the importance of considering multiple physicochemical processes (precipitation/dissolution, organic complexation, adsorption/desorption and ionic association dissociation), bio-uptake and components (Fe, Ni, Co, Ca, Mg, Se, Mo) as responsible for TE dynamics in AD. To the best of our knowledge, there is no AD mathematical model which considers TE dynamics and attempts to quantify the complete profile of TE speciation across different phases. Therefore, in order to accurately predict the dynamics of TEs in AD, inclusion of additional processes and variables is required to the ADM1 model.

The main objective of this work is to develop a modelling framework to: (1) mechanistically describe all the main biochemical and physicochemical processes involved in an AD system, including adsorption of TEs (on inerts, biomass and precipitates), EDTA/VFA complexation, precipitation/dissolution and uptake of TEs; (2) investigate the speciation of TEs in AD, in particular the partitioning of TEs among adsorbed species, precipitates and dissolved complexes in relation to the change in initial concentration of complex organic matter, TEs, calcium, magnesium and EDTA. To this aim, new TE-adsorption reactions have been defined and added to the ADM1-based model developed<sup>37</sup>. Hydrolysis, acidogenesis, acetogenesis, methanogenesis and acid-base equilibrium have been combined to TE-adsorption, complexation and precipitation/dissolution and the effects on total methane production have been evaluated. Further, experimental data from a real case study have been compared with the model predictions in continuous mode of operation.

## Results

The implementation of the developed mathematical model has involved: (i) the update of the biochemical framework to include the effect of TEs on microbial kinetics as well as the release of sulfur, phosphorous and TEs during the disintegration step; (ii) the introduction of new acid-base reactions; (iii) inclusion of new liquid-liquid processes (TE-EDTA and TE-VFA complexation); (iv) the inclusion of solid-liquid processes (precipitation and dissolution); (v) incorporation of dynamic state variables and processes for sorption/desorption processes; (vi) the collection of the parameters values governing the sorption processes. In particular, these values have been collected from suitable literature and scientific reports<sup>42</sup>.

The capability of the model to cover a wide range of input variables and parameters is crucial to model applicability. In this regard, five different simulation scenarios (Table 1) have been considered to check the model prediction capability. *Scenario 1* analysed the impact of an increase in adsorbent concentration (complex organic matter) when the amount of TEs in the system is constant. *Scenario 2* evaluated the change in speciation of TEs when the initial concentration of a single TE is changed. Here the model simulations were performed for Fe, Ni and Co. *Scenario 3* examined any possible competition among carbonate and sulfide systems to capture mobile TEs onto available surfaces. *Scenario 4* was used to assess the response of the sorption reactions to the presence of a synthetic chelating agent. Lastly, *Scenario 5* showcased the variation in the model prediction due to change in binding site density for biomass, inerts and precipitates. In addition to *in-silico* experiments, a continuous mode of operation was implemented to compare model behaviour with a real case study demonstrating the effect of TE starvation on methane production.

**Scenario 1: Effect of adsorbents on TE dynamics in anaerobic digestion.** *Scenario 1* assesses the change in TE sorption with increase in $X_c$ amount in the reactor with $B_δ = 2 \times 10^{-3}$, $I_δ = 2 \times 10^{-3}$ and $P_δ = 2 \times 10^{-5}$ moles/g of dry weight. $B_δ$ is the biomass binding site density, $I_δ$ is the inert binding site density and $P_δ$ is the precipitate binding site density. The initial TEs concentration and other initial concentrations were kept constant (Tables 1, 2). What stands out in this scenario is the increase in adsorption of TEs to inert with increase in initial $X_c$ (Fig. 1). This is more pronounced in Fe fractions where $I-Fe_{(s)}$ increases from around 45% of Fe to about 70% of Fe for corresponding initial $X_c$. Ni and Co also share similar behaviour with inert. Consequently,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** TE speciation and methane production in an anaerobic batch reactor, under different initial complex organic matter concentrations.</figcaption></figure>

the fraction for TE-EDTA decreased with increase in initial $X_c$. Similar is the case for biomass fraction, $B-TE_{(s)}$. Further, with the increase in initial $X_c$, the sulfide precipitate decreased. Sulfide fraction for Fe was higher than Ni and Co for corresponding initial $X_c$. The majority of Ni and Co (< 70%) was found associated to EDTA. For Ni and Co, the TE-EDTA fraction remained more or less constant across all initial $X_c$ concentrations (1–5 gCOD of composite organic matter). It is interesting to note (Tables 1, 2) that for all model runs, the cumulative methane production (Fig. 1) increased. This may be due to increased amount of organic matter and an increased TEs release into the system which is proportional to $X_c$.

**Scenario 2: Effect of initial Fe, Ni and Co concentration on TE dynamics in anaerobic digestion.** Scenario 2 studies the change in TE dynamics due to change in initial concentrations of TEs. Model simulation for 8 different initial concentrations for each TE (Fe, Ni and Co) was carried out (Fig. 2a–c). The initial concentration for the TE under study was varied while the other two TEs concentrations were set at an optimal value (Tables 1, 2).

*Change in initial $S_{Fe^{2+}}$ concentration.* The amount of Fe adsorbed on inert matter is found to be higher at both lower and higher initial concentration of $S_{Fe^{2+}}$(Fig. 2a). In RUN 1 through RUN 8, inert fraction of TEs is higher for Fe. Ni and Co have comparatively lower inert fraction, which increases rapidly in comparison to Fe inert fraction. Change in initial concentration of $S_{Fe^{2+}}$ resulted in higher sulfide fraction for Fe and Co as compared to Ni. Ni and Co were found more with EDTA complex fraction. Additionally, with increase in initial $S_{Fe^{2+}}$ concentration the carbonate fraction increases. The carbonate fraction for RUN 8 (with highest initial $S_{Fe^{2+}}$ concentration, $3 \times 10^{-3}$ M) is the highest among all the runs irrespective of the TE fraction. It is also interesting to note that complexation fraction is only observed in case of Ni and Co which decreases with increase in $S_{Fe^{2+}}$ concentration. The free TE fraction for Ni and Co increases slowly with increase in $S_{Fe^{2+}}$ initial concentration but is found to be constant for Fe. The cumulative $CH_4$ production profile differs only for RUN 8 which corresponds to the highest $S_{Fe^{2+}}$ concentration.

*Change in initial $S_{Ni^{2+}}$ concentration.* TE-inert fraction is predominant in Fe for all initial concentrations of $S_{Ni^{2+}}$ (Fig. 2b). However, for Ni and Co fractions, EDTA complex fraction is higher at lower initial concentration of $S_{Ni^{2+}}$. With increase in initial $S_{Ni^{2+}}$ concentration, the inert fraction increases for Ni and Co. Sulfide fraction was more prominent for Fe and Co. In case of Fe speciation, the fraction for sulfide and inert remains more or less constant for all the values of $S_{Ni^{2+}}$ initial concentration. A thorough analysis of the fractionation pattern shows that with increase in initial $S_{Ni^{2+}}$ concentration, the carbonate fraction does not increase significantly for Fe and Ni. It is also interesting to note here that complexation fraction decreases with increase in $S_{Ni^{2+}}$ initial concentration. For Fe, the complexation fraction decreases with higher $S_{Ni^{2+}}$ initial concentration. The free TE fraction changes with increase in $S_{Ni^{2+}}$ initial concentration for Co and for Ni. In case of Fe, the free TE remained more or less constant throughout all the simulations. RUN 5 showed the highest cumulative methane production. Subsequently, the methane production decreased with further increase in initial $S_{Ni^{2+}}$ concentration. The cumulative methane production for RUN 1, with $5 \times 10^{-7}$ M of Ni, is around $2.3 \times 10^{-2}$ M. For RUN 8 methane production is the lowest among all simulation runs.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) TE speciation and methane production in an anaerobic batch reactor, under different initial Fe concentrations. (b) TE speciation and methane production in an anaerobic batch reactor, under different initial Ni concentrations. (c) TE speciation and methane production in an anaerobic batch reactor, under different initial Co concentrations.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (continued)</figcaption></figure>

*Change in initial S<sub>Co2+</sub> concentration.* The trend of fractionation among Fe, Ni and Co is similar to the previous simulation sets for change in initial S<sub>Fe2+</sub> and S<sub>Ni2+</sub> concentrations (Fig. 2c). Higher amount of Fe was found adsorbed on inert matter at all initial concentration of S<sub>Co2+</sub>. Comparatively, higher amount of sulfide was observed in Fe and Co fractions, and EDTA complexation fraction was higher for Ni. A very minute fraction of complexation is observed in all model runs for Fe. The fractionation of Fe shows little or no change with change in initial S<sub>Co2+</sub> concentration. The carbonate fraction does not increase significantly for all the TEs. The free TE fraction does not change significantly for any of TEs. The cumulative methane production did not change drastically with change in initial S<sub>Co2+</sub> concentration. The RUNs (1–8) resulted in cumulative methane production close to around 2.25 × 10<sup>-2</sup> M. However, there was a variation in the rate of methane production.

**Scenario 3: Effect of initial S<sub>Ca2+</sub> and S<sub>Mg2+</sub> concentration on TE dynamics in anaerobic digestion.** Scenario 3 investigates the effect of change in initial S<sub>Ca2+</sub> and S<sub>Mg2+</sub> concentrations on speciation of TEs and cumulative methane production (Fig. 3). With increase in S<sub>Ca2+</sub> and S<sub>Mg2+</sub> initial concentrations, inert and biomass fractions decreased for Fe, Ni and Co. For example, the amount of Fe bound to inert matter at lower initial concentration of S<sub>Ca2+</sub> and S<sub>Mg2+</sub> is higher in RUN 1 as compared to RUN 8. Conversely, amount of carbonate fraction for all the TEs increased with increase in S<sub>Ca2+</sub> and S<sub>Mg2+</sub> initial concentration. Free TE fraction was only noticed for Fe. For Ni and Co, the free TE fraction was close to zero and only appeared in the runs with lower initial concentration of S<sub>Ca2+</sub> and S<sub>Mg2+</sub>. Sulfide fraction is principally found with Fe and increases for all TEs with increase in initial S<sub>Ca2+</sub> and S<sub>Mg2+</sub> concentration. EDTA complex fraction was found across all the simulations for Ni and Co. Fe was found associated more with inert fraction as well as biomass fraction with lower S<sub>Ca2+</sub> and S<sub>Mg2+</sub> concentration. With the increase in S<sub>Ca2+</sub> and S<sub>Mg2+</sub> initial concentration rate of cumulative CH<sub>4</sub> production decreased, while the total amount of CH<sub>4</sub> produced remained more or less constant for RUN 1–7 with a significant decrease for RUN8 which had highest amount of initial S<sub>Ca2+</sub> and S<sub>Mg2+</sub> concentration.

**Scenario 4: Effect of initial EDTA concentration on TE dynamics in anaerobic digestion.** Scenario 4 explored the effect of initial EDTA concentrations on the sorption of TEs and cumulative CH₄ production with a set of 10 simulations (Fig. 4). The initial concentrations of Fe, Ni and Co are set at optimum values for all simulations (Tables 1, 2). The initial concentration of EDTA was varied from 1 × 10⁻⁵ M to 1 × 10⁻⁴ M (Tables 1, 2). Fe was found predominantly in inert fraction with lower initial concentration of EDTA. For example, the Fe-inert fraction is higher at RUN 1 as compared to RUN 10. As the initial concentration of EDTA increased, the Fe was observed more in the complexed fraction. However, Ni and Co were primarily found with complex fraction, even with lower initial EDTA concentration. It is interesting to note that, Fe is also found with sulfide fraction, whereas, Ni and Co is found with sulfide and inert fraction only at lower initial concentration of EDTA. The fraction of Fe bound to inert is higher as compared to the sulfide and complexed fraction for lower initial EDTA concentration. The free form of TEs decreased with increase in EDTA initial concentration. S<sub>Fe2+</sub> decreased with increase in EDTA initial concentration and a small fraction (~5% Fe) was found in the free form in RUN 10. However, Ni and Co free fractions were not observed in the model simulations other than RUN 1. Addition of higher amount of initial EDTA into the system did not significantly affect cumulative CH₄ pro-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** TE speciation and methane production in an anaerobic batch reactor, under different initial Ca-Mg concentrations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** TE speciation and methane production in an anaerobic batch reactor, under different initial EDTA concentrations.</figcaption></figure>

duction. Nevertheless, the rate of methane production increased with higher initial EDTA concentration. The amount of methane produced at RUN 1 is close to $2.45 \times 10^{-2}$ M and the amount of methane produced during RUN 10 is approximately $2.5 \times 10^{-2}$ M.

**Scenario 5: Effect of binding site density on TE dynamics in anaerobic digestion.** *Change in biomass binding site density* $B_δ$. Nine simulations were carried out with different binding site density ($B_δ$) for biomass (Fig. 5a, Tables 1, 2). The different binding site density was chosen while keeping in mind the binding site density reported in literature, which is $2 \times 10^{-3}$ mol/g of biomass. With decrease in $B_δ$ values, the fraction of Fe bound to biomass decreases (RUN 1 through 9, Fig. 5a). A similar trend has been observed for Co. This is also true for FeS and CoS. Change in $B_δ$ values could affect the speciation of TEs but the distribution of TEs across
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) TE speciation and methane production in an anaerobic batch reactor under different biomass binding site density. (b) TE speciation and methane production in an anaerobic batch reactor under different inert binding site density. (c) TE speciation and methane production in an anaerobic batch reactor, under different precipitate binding site density.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** (continued)</figcaption></figure>

different fractions remained more or less constant for RUN 5 through RUN 9. Maximum changes in free TE fraction were observed between RUN 1 and 5. With further decrease in $B_δ$ values the speciation pattern did not change. Similarly, cumulative $CH_4$ decreased till RUN 5 after which it remained constant.

*Change in inert binding site density I*<sub>δ</sub>. Figure 5b depicts 9 simulations carried out with 9 different $I_δ$ values (Tables 1, 2). Change in $I_δ$ altered the inert fraction for all model RUNs. The inert fraction for Fe decreased as $I_δ$ decreased from RUN 1 to RUN 9. Whereas, the increase in binding site density for inert resulted in an increase in sulfide fraction. For Co and Ni similar speciation was observed but with two prominent characteristics. One, the free TE fraction decreased with decrease in inert binding site density for Ni and Co whereas in case of Fe, the free TE fraction increased up to RUN 5 and then decreased in RUN 9. Two, the EDTA complex fractions for Ni and Co were comparatively higher as compared to EDTA complexation fraction of Fe. The cumulative $CH_4$ production was also examined (Fig. 5b) with change in $I_δ$. RUN 9, with lowest inert binding site density, resulted in highest cumulative $CH_4$ production.

*Change in precipitate binding site density P*<sub>δ</sub>. Figure 5c presents 8 simulations carried out with 8 different binding site densities for precipitates. In this model only FeS has been considered to harbour TE binding sites. The binding site density of FeS for Ni and Co has been assumed in this study (see "Change in precipitate binding site density P*δ*"). The changes in binding site density of mineral precipitate did not result in change of TE speciation. The speciation for Fe, Ni and Co remains intact all throughout the RUNs from 1 to 8. Same is the case with cumulative methane production.

**Model comparison with a real case study.** The model output has been qualitatively compared with the experimental results of a highly cited work<sup>43</sup>. The study assessed the impact of TE supplementation on single stage continuous mesophilic AD of food waste operated at bench scale. The substrate was a grounded mixture of fruits, vegetable, meat, fried foods and fat trimmings. A constant organic loading rate of 1.45 g/VSl/day was maintained for the reactors. The reactors were operated for 120 days at multiple HRT (25, 50, 100 and 180 (days)) with two reactors for each HRT. Each pair of such reactors was run with the same feedstock characteristics, but one of these reactors received TE supplementation periodically while the other one did not. The TEs have been dosed in a semi-continuous operation mode. The authors highlighted the failure of the reactors, after around 100 days, with no TE supplementation, whereas the reactors with TE supplementation exhibited a stable biogas production.

The starvation reactor with 50 days HRT was chosen for comparing the model output in terms of biogas production. The substrate characteristics and the operation parameters have been reported in Table 3. The model parameters, coefficients and stoichiometric values, have been selected based on the ADM1. Any change in the benchmarked model parameters have been explicitly mentioned in the text. The coefficients of the dose response function for a specific TE have been set at suitable values to capture starvation effect of TEs properly.

The comparison of experimental data and model predictions for methane production has been summarized in Fig. 6a. According to the experimental results, TEs (Fe, Ni and Co) were the limiting factor which caused reactor failure after 100 days of operation. The methane production fell from apparently steady methane production
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** (a) Comparison of experimental data and model predictions for methane production of an AD reactor operated continuously with a 50 days HRT and food waste as feedstock. (b) Simulated dynamics of free TEs, precipitates and speciation of sulfide during model comparison with real case data.</figcaption></figure>

of 0.43 L/gVS to 0 in the last 20 days of operation (Fig. 6a). The model has satisfactorily predicted the methane production during the AD of food waste. In particular, the model was able to reproduce the starvation effect of TEs on methane production in a continuous AD system. Furthermore, the free TEs along with speciation have been plotted in Fig. 6b, where TE depletion has been reported along with corresponding dynamics of sulfide precipitate in the system. Depletion of Ni, Co and Fe from the reactor have been the limiting factor which resulted in a reactor failure after around 100 days, which is the major finding of the experimental result$^{43}$.

Discussion

In an AD system, depending upon the reactor conditions (i.e. pH, OLR, redox state and hydraulic residence time), TEs are found in different chemical forms or species. The fate of TEs depend primarily on various biogeo-
chemical processes occurring simultaneously at different rates. Precipitation, dissolution, organic complexation,
adsorption, biotransformation of TE containing intermediates and bio-uptake are such determining processes.
The resultant amount of free TE which eventually remains in the liquid phase depends on the amount of TEs
involved in these mentioned biogeochemical fluxes. Precipitation/dissolution is a major flux which results in
formation of sulfide, carbonate and phosphate minerals, and consequently scavenges the majority of TEs from
the system. Similarly, organic complexation utilizes a portion of free TEs. TEs end up complexed with natural
or synthetic organic ligands present in AD systems. A fraction of TEs is up-taken by biological cells and hence
initiates other cascading processes (such as gene regulation, storage, excretion)$^{19,44}$. However, some amount of
TEs is supplied into the system by the biodegradation (in our case disintegration of complex organic matter) of
TE containing intermediates. Adsorption/desorption of TE is an important and major physicochemical process,
particularly when a large amount of surface area (biomass, inorganic precipitate and inert material) is avail-
able in the AD systems. The proposed model under discussion includes all the significant processes (including
adsorption) involved in the speciation of TEs in AD. Apart from the standard biodegradation reactions (as in
ADM1) the model includes 42 adsorption/desorption reactions, 13 precipitation/dissolution reactions and 15
complexation reactions. The adsorption of TEs has been considered on the biomass, precipitate and inert matter.
The proposed model represents a general framework wherein the user/modeller can add and remove relevant
inorganic components and species of interest as well as corresponding processes, to dynamically predict the TEs
speciation. Such a customizable general model can be used as a precursor to build various model-based applica-
tions that aims at quantifying the fate of TEs in anaerobic digestion environments.

The proportional variations and dynamics of chemical fractions for adsorbed TEs under different initial concentration of: (i) adsorbents (Fig. 1), (ii) TE (Fe, Ni and Co, Fig. 2), Ca-Mg (Fig. 3), EDTA (Fig. 4) and the effect of change in binding site density (Fig. 5) have been studied. Throughout all the simulation scenarios, the adsorption fraction for inerts is higher than the adsorption fraction for biomass and precipitates. A possible explanation for this observation is that amount of inert material formed by the end of AD process is always higher as compared to the amount of biomass and mineral precipitate. Secondly, the % fraction is calculated at the end of the 100 days when the amount of biomass is already at lower limits and the inert materials reach maximum values. Thirdly, the binding site density may also influence the adsorption fraction.

Each simulation scenario is discussed individually below. Due to scarcity of literature we have tried to offer
possible explanations based on theoretical understandings, first principles and our experience.

Change in initial concentration of Xc, as reported in Fig. 1, affects the TE adsorption fraction. The increase in TE adsorption fraction is more visible for inerts. There is no significant increase in adsorption to biomass and precipitate. This is primarily because with increase in Xc the amount of inerts increases rapidly and accumulates in the system more efficiently as compared to biomass and precipitates. With increase in adsorbent concentration in the system the amount of adsorbed TEs increases (particularly with inert matter) for Fe, Ni and Co. This is because increase in adsorbent amount provides greater number of surface binding sites and hence the increase in adsorption. The decrease in free TEs can be attributed to the uptake of TEs by microbial communities as well as increase in adsorption flux.

The main objective of the scenarios with change in initial Fe (Fig. 2a), Ni (Fig. 2b) and Co (Fig. 2c) was to predict the adsorption behaviour when limited and over supplied (in comparison to optimal concentration values) amount of TEs are present in the AD system. The simulation results indicate that the initial concentrations of Fe, Ni and Co have distinct patterns of effect on adsorption. Adsorption pattern as seen in case of different initial Fe concentrations (Fig. 2a) implies a decrease in adsorption fraction (for Fe and Co) with increase in TE dose. This may be explained by the fact that an increase in TE concentrations implies that TEs end up in sulfide fraction due to a suitable pH regime. Also, the difference in cumulative methane production for lower Fe concentration (RUN 1) and higher Fe concentration (RUN 8) is small because addition of more Fe does not affect the system. There may be two possible explanations: (a) as per our experience and literature anaerobic digestion system is more sensitive to change in Ni concentration; and (b) TE-EDTA complexes play a major role in stimulating methane production. However, this is not the case with change in initial Ni concentration (Fig. 2b) where a small change in initial concentration significantly affects the cumulative methane production. This may be because there is already enough free Ni in the system. Addition of a small amount of Ni into the system greatly affects the system both in terms of methane production and speciation of TEs. The increase in sulfide fraction and decrease in inert fraction implies an underlying mass transfer from inert fractions to the sulfide fraction. It can also be postulated that first the TEs get adsorbed on the surfaces and slowly get precipitated with increase in Ni dosing. Unlike Fe and Ni, change in initial Co concentration did not translate into any significant change in the TE speciation pattern (Fig. 2c). This may be attributed to the low concentration and high solubility of Co in the reactor system which can alter the reactor equilibrium significantly. A similar observation was made in a study²⁰ which concluded an increased solubility of Co in presence of sulfide precipitates along with Fe and Ni. Interestingly, the model is able to capture small changes in methane production due to change in Co concentrations. This is an important aspect of the model, particularly when sensitivity of methane production to small concentration of TEs is considered.

Ca and Mg have been reported to exhibit competitive effect with TEs in terms of occupying binding sites. This is because of the similar charge of Ca and Mg as divalent TE ions. The decrease in methane production (Fig. 3) due to increase in Ca-Mg addition (Tables 1, 2), can be ascribed to the change in pH of the system and the availability of TEs. Moreover, the decrease in adsorption fraction with increase in Ca-Mg is due to increase in sulfide precipitate. At higher concentration of Ca-Mg, TEs precipitate with carbonates rather than adsorb on inert. This can be explained by the increase in saturation level of calcium-carbonate-magnesium systems due to increase in Ca-Mg concentration. This exercise confirms that with increase in TE dose the cumulative methane production first increases then decreases when the dose amount exceeds the optimal concentration limit. However, with increase in Ca-Mg dosage there is only decrease in cumulative methane production. This can partly be explained by the decrease in availability of free TEs at the higher concentration of Ca-Mg dosage.

Change in initial EDTA concentration, as reported in Fig. 4, affects the TE adsorption fraction. The decrease in TE adsorption fraction with addition of EDTA confirms its chelating ability. The inert fraction decreased along with decrease in sulfide precipitate fraction. The decrease in sulfide precipitate is in accordance to our experience³⁷. The modelling framework presented here can be modified to take into account the chemistry of green chelating agents such as EDDS which has recently been investigated for its application and fate in anaerobic systems⁴¹,⁴⁵,⁴⁶.

Model application is not straightforward in the case of AD systems. Identifying correct inlet stream and operational parameters is crucial to get accurate model predictions. Consequently, model simulation depends on availability of data sets with all the necessary input characteristics, in particular the characterization of organic
waste (particulate carbohydrate, proteins and lipids) in terms of COD and with regard to TEs. We selected a particular data set from published reports on TE supplementation in AD of food waste. The data set was selected by keeping in mind the specific effect of TEs (here starvation) it demonstrates in the continuous reactor operation.

Model calibration and validation is more complicated when physicochemical processes such as precipitation, complexation and adsorption are considered in addition to biodegradation processes. This is because the range of TEs concentration, in terms of free TEs, required for optimum microbial activity is rarely reported in conjugation with methane production and other major intermediates of AD process. Moreover, in most studies organic matter is characterized in terms of volatile solids content and the relationship between volatile solids content and COD varies for each specific substrate. Also, the characterization of organic waste in terms of content of sulfur, phosphorous and TEs is difficult to find for a particular study simultaneously. The primary reason for such a discrepancy is due to the variety of substrates used and the diversity of research objectives. Nevertheless, with the increasing amount of data in recent years dedicated to microbial community analysis in AD, the effect of TEs on methanogenesis and experimental speciation analysis, it would be relatively more convenient to calibrate and validate TE dynamics and speciation in AD system in the near future.

The recruitment of a model-based control strategy of TEs addition to AD system may require additional model complexity in terms of significant processes and components. Likewise, considering EPS (extracellular polymeric substances) as a storage house for complexation of TEs, which in turn affects the binding site density and hence adsorption behaviour, it is necessary to incorporate EPS and SMP (soluble microbial products) in order to achieve a more accurate mathematical model. Nevertheless, the goal of this study is to build a model framework based on ADM1 which can simulate the effect of TEs dynamics on AD. This model can be used to formulate a TE dosing scheme or a reduced version of this can be used as a model predictive control of TE dosing in AD. At this stage, with unavailability of proper experimental data to calibrate and validate such a model, it is necessary to individuate and define other biochemical and physicochemical processes (EPS and redox state) which can affect the speciation of TEs in anaerobic digesters. The usefulness of such a complex dynamic model based on ADM1 to predict the speciation of TEs in anaerobic digestion system may be questioned. Nevertheless, the need to implement major processes affecting of TE dynamics in a consolidated framework is crucial for the preceding steps of model development. Therefore, some features of the model are under development and some limitations and gap areas in this study are summarized below:

1. Amino acids and proteins can be considered as a source for TEs.

2. Production of soluble microbial product as a part of biochemical processes can be included.

3. Incorporating a separate module to define and predict binding site density by considering surface complexation models.

4. Information on microbial community composition and abundance in presence of TEs can be incorporated.

Future developments include a global sensitivity analysis to individuate the most significant parameters; model calibration and validation by using ad hoc experimental data.

## Methods: model development and simulation

**Mathematical model structure.** The proposed model simulates the dynamics of TEs during AD (Fig. 7). The 3-phase (solid-liquid-gas) AD model, which is based on mass conservation principle, considers a full kinetic framework to incorporate processes responsible for speciation of TEs. Speciation reactions occurring in the liquid-solid phase are central to the TE adsorption and they can be classified in three groups: adsorption on biomass, adsorption on inerts and adsorption on TE precipitates (i.e. FeS). Apart from adsorption reactions, protonation/deprotonation, TE precipitation/dissolution and TE-complexation with VFAs and EDTA have been considered³⁷. Biochemical degradation of organic complex matter in the liquid phase has been considered to produce a gas mixture of H₂, CH₄, CO₂ and H₂S.

Overall, the proposed model tracks the dynamics of the state variables which represent the components of the proposed AD model. For each component considered in the liquid, gas and solid phase the general mass balance can be written as:

$$ \frac{dV_{\text{liq}}S_i}{dt} = q_{\text{in}}S_{\text{in},i} - q_{\text{out}}S_i + V_{\text{liq}} \left( \sum_{j=1}^{m_5} \delta_{i,j} \rho_{A,j}(t, \mathbf{S}) - \rho_{T,i}(t, \mathbf{S}, \mathbf{S}_{\text{gas}}) \right) \\ + V_{\text{liq}} \left( \sum_{j=1}^{m_1} \alpha_{i,j} \rho_{\text{bio},j}(t, \mathbf{S}, \mathbf{X}) + \sum_{j=1}^{m_2} \beta_{i,j} \rho_{\text{cmplx},j}(t, \mathbf{S}) \right) \\ + V_{\text{liq}} \left( \sum_{j=1}^{m_3} \gamma_{i,j} \rho_{\text{prec},j}(t, \mathbf{S}, \mathbf{X_p}) + \sum_{j=1}^{m_3} \gamma_{i,j} \rho_{\text{dissol},j}(t, \mathbf{S}, \mathbf{X_p}) + \sum_{j=1}^{m_4} \varepsilon_{i,j} \rho_{s,j}(t, \mathbf{S}, \mathbf{X}_{\text{ads}}, \bar{\mathbf{X}}_{\text{ads}}) \right), \quad (1) \\ i = 1, \dots, n_1, \quad t > 0, $$
<figure><img src="image_12.png" /><figcaption>Figure 7. The TE-ADM1 model. Created in Lucidchart (www.lucidchart.com).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8. A simplified presentation of the sorption/desorption model. rho1 represents the disintegration and hydrolysis step for degradation of Xc to simpler monomers such as Ssu, Saa, Sfa, XI and S1. rho2 is the set of all subsequent processes in AD which lead to biomass growth and subsequent degradation of soluble metabolite into CH4. rho3/rho4 is the new adorption desorption process in ADM1 where FB_Xsu and OB_Xsu. represent the free and binding sites concentration on biomass Xsu. rho5 is the biomass decay and rho6 is the release of TEs concomitant to biomass decay. Created in Lucidchart (www.lucidchart.com).</figcaption></figure>

$$ \frac{dV_{\text{liq}}X_i}{dt} = q_{\text{in}}X_{\text{in},i} - q_{\text{out}}X_i + V_{\text{liq}} \left( \sum_{j=1}^{m_1} \alpha_{i,j} \rho_{\text{bio},j}(t, \mathbf{S}, \mathbf{X}) \right), \quad (2) $$

$$ \frac{dS_{\text{gas},i}}{dt} = - \frac{S_{\text{gas},i}q_{\text{gas}}}{V_{\text{gas}}} + \frac{V_{\text{liq}}}{V_{\text{gas}}} \rho_{\text{T},i}(t, \mathbf{S}, \mathbf{S}_{\text{gas}}), \quad (3) $$
<table><caption>**Table 1.** Variable initial concentration of complex organic matter, TEs, EDTA, Ca and Mg used in the model simulations.</caption><thead><tr><th rowspan="2">RUN</th><th rowspan="2">Scenario 1<br>Variable COM<br>(gCOD)</th><th colspan="3">Scenario 2<br>Variable TEs (M)</th><th rowspan="2">Scenario 3<br>Variable Ca-<br>Mg (M)</th><th rowspan="2">Scenario 4<br>Variable EDTA (M)</th><th colspan="3">Scenario 5<br>Variable Binding density</th></tr><tr><th>Fe</th><th>Ni</th><th>Co</th><th>Biomass</th><th>Inert</th><th>Precipitate</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>5.0 × 10<sup>-6</sup></td><td>5.0 × 10<sup>-7</sup></td><td>5.0 × 10<sup>-9</sup></td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>1</sup></td><td>2.0 × 10<sup>1</sup></td><td>2.0 × 10<sup>-1</sup></td></tr><tr><td>2</td><td>2</td><td>9.0 × 10<sup>-6</sup></td><td>9.0 × 10<sup>-7</sup></td><td>9.0 × 10<sup>-9</sup></td><td>1.5 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-5</sup></td><td>2.0</td><td>2.0</td><td>2.0 × 10<sup>-2</sup></td></tr><tr><td>3</td><td>3</td><td>3.0 × 10<sup>-5</sup></td><td>3.0 × 10<sup>-6</sup></td><td>3.0 × 10<sup>-8</sup></td><td>2.0 × 10<sup>-3</sup></td><td>3.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-1</sup></td><td>2.0 × 10<sup>-1</sup></td><td>2.0 × 10<sup>-3</sup></td></tr><tr><td>4</td><td>4</td><td>7.0 × 10<sup>-5</sup></td><td>7.0 × 10<sup>-6</sup></td><td>7.0 × 10<sup>-8</sup></td><td>2.5 × 10<sup>-3</sup></td><td>4.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-4</sup></td></tr><tr><td>5</td><td>5</td><td>1.0 × 10<sup>-4</sup></td><td>1.0 × 10<sup>-5</sup></td><td>1.0 × 10<sup>-7</sup></td><td>3.0 × 10<sup>-3</sup></td><td>5.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-5</sup></td></tr><tr><td>6</td><td></td><td>5.0 × 10<sup>-4</sup></td><td>5.0 × 10<sup>-5</sup></td><td>5.0 × 10<sup>-7</sup></td><td>3.5 × 10<sup>-3</sup></td><td>6.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-6</sup></td></tr><tr><td>7</td><td></td><td>9.0 × 10<sup>-4</sup></td><td>9.0 × 10<sup>-5</sup></td><td>9.0 × 10<sup>-7</sup></td><td>4.0 × 10<sup>-3</sup></td><td>7.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-7</sup></td></tr><tr><td>8</td><td></td><td>3.0 × 10<sup>-3</sup></td><td>3.0 × 10<sup>-4</sup></td><td>3.0 × 10<sup>-6</sup></td><td>4.5 × 10<sup>-3</sup></td><td>8.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-6</sup></td><td>2.0 × 10<sup>-6</sup></td><td>2.0 × 10<sup>-8</sup></td></tr><tr><td>9</td><td></td><td></td><td></td><td></td><td></td><td>9.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-7</sup></td><td>2.0 × 10<sup>-7</sup></td><td></td></tr><tr><td>10</td><td></td><td></td><td></td><td></td><td></td><td>1.0 × 10<sup>-4</sup></td><td></td><td></td><td></td></tr></tbody></table>

$$ \frac{dV_{\text{liq}}X_{\text{p,i}}}{dt} = q_{\text{in}}X_{\text{p}}^{\text{in}} - q_{\text{out}}X_{\text{p,i}} + V_{\text{liq}}(\bar{\alpha}_{i}\rho_{\text{prec,i}}(t, \mathbf{S}, \mathbf{X}_{\text{p}}) + \bar{\beta}_{i}\rho_{\text{dissol,i}}(t, \mathbf{S}, \mathbf{X}_{\text{p}})), \quad (4) $$
$$ i = 1, \dots, m_3, \quad t > 0, $$

where: $n_1$ is the number of soluble components; $n_2 - n_1$ is the number of particulate components, including complex organic matter and inert; $n_3 = n_2 - n_1 + m_3$ is the number of adsorbing components, including microbial species, inerts and precipitates; $m_1$ is the number of biochemical processes considered; $m_2$ is the number of complexation processes considered; $m_3$ is the number of precipitation/dissolution processes considered; $m_4$ is the number of adsorption/desorption processes considered; $m_5$ is the number of acid-base reactions considered; $\alpha_{i,j}$ is the stoichiometric coefficient of the $i$th species on biochemical process $j$; $\beta_{i,j}$ is the stoichiometric coefficient of the $i$th species on complexation process $j$; $\gamma_{i,j}$ is the stoichiometric coefficient of the $i$th species on precipitation/dissolution process $j$; $\delta_{i,j}$ is the stoichiometric coefficient of the $i$th species on acid-base reaction $j$; $\epsilon_{i,j}$ is the stoichiometric coefficient of the $i$th species on adsorption/desorption reaction $j$; $\bar{\alpha}_i$ and $\bar{\beta}_i$ are the stoichiometric coefficients for the precipitate formation and dissolution; $S_i$ is the $i$th soluble component, $S = (S_1, \dots, S_{n_1})$; $X_i$ is the $i$th particulate component, $X = (X_{n_1+1}, \dots, X_{n_2})$; $X_{p,i}$ is the $i$th precipitate, $X_p = (X_{p,1}, \dots, X_{p,m_3})$; $X_{ads,i}$ is the concentration of free binding sites for a specific adsorbent component, $X_{ads} = (X_{ads,1}, \dots, X_{ads,n_3})$; $\bar{X}_{ads,i}$ is the concentration of occupied binding sites for a specific adsorbent component, $\bar{X}_{ads} = (\bar{X}_{ads,1}, \dots, \bar{X}_{ads,n_3})$; $S_{gas,i}$ is the $i$th soluble gas component, $S_{gas} = (S_{gas,1}, \dots, S_{gas,n_1})$; $\rho_{A,j}(t, S)$ is the process rate for the $j$th acid base reaction; $\rho_{T,i}(t, S, S_{gas})$ is the gas transfer process rate for the $i$th soluble component; $\rho_{bio,j}(t, S, X)$ is the rate for the $j$th biochemical process; $\rho_{complx,j}(t, S, X)$ is the rate of the $j$th complexation process; $\rho_{prec,j}(t, S, X_p)$ is the rate for the $j$th precipitation process; $\rho_{dissol,j}(t, S, X_{ads}, \bar{X}_{ads})$ is the rate for the $j$th dissolution process; $\rho_{s,j}(t, S, X_{ads}, \bar{X}_{ads})$ is the rate for the $j$th sorption/desorption process; $V_{\text{liq}}$ and $V_{\text{gas}}$ are the liquid and gas volume of the reactor; $q_{\text{gas}}$ is the gas flow rate.

**Biochemical fate of TE during AD.** During AD, microorganisms are able to uptake a limited amount of TEs available in the bulk phase. However, the original ADM1 does not consider TEs as state variables by altogether neglecting any inorganic component other than ammonium. Starting with complex organic matter as the first reactant/component, AD is considered to have 5 stages of biochemical reactions: disintegration and hydrolysis which constitute extracellular processes; acidogenesis, acetogenesis and methanogenesis as intracellular processes. The organic complex matter is assumed to be a critical source of simpler compounds, as it eventually disintegrates to form carbohydrates, proteins and lipids. The model proposed in this paper considers the release of TEs and inorganic components during the disintegration process. Phosphorus, sulfur and alkali earth elements have been set as a product of the disintegration step. Phosphorous has been supposed to be released as $HPO_4^{2-}$, which represents the most abundant form in the pH range of 6–14. Likewise, the released sulfur has been considered in the form of $HS^-$ which has been depicted as the most abundant form in the pH range of 6–12. Moreover, external dosing of TEs is also included in the model as additional TEs source. TEs bio-uptake has been linked to acetate and hydrogen uptake. The biochemical stoichiometry has been updated accordingly.

**Acid-base processes affecting TEs geochemistry.** In addition to the 6 acid base reactions of the original ADM1, the model includes sulfate acid base system, phosphate acid base system and carbonate acid base system to define the effect of pH and other biochemical processes on the geochemistry of TEs. The new inorganic components include bicarbonate/carbonate, phosphoric acid/dihydrogen phosphate, dihydrogen phosphoric acid/hydrogen phosphate, hydrogenphosphoric acid/phosphate, hydrogen sulfide/bisulfide and bisulfide/sulfide acid pairs. Apart from inorganic components, acid base reactions for synthetic chelators (here EDTA) have also been added as they are involved in TEs complexation reactions. The corresponding thermodynamic values and rate coefficients for the acid base processes have been sourced from literature. Charge balance has been modified in order to consider the effect of additional components which may affect pH prediction.
<table><caption>**Table 2.** Initial concentration of various dynamic state variables considered in the model for the batch studies.

ªRefer Table 1.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Variable</th><th rowspan="2">Scenario 1 COM</th><th colspan="3">Scenario 2 Variable TEs</th><th rowspan="2">Scenario 3 Ca-Mg</th><th rowspan="2">Scenario 4 EDTA</th><th colspan="3">Scenario 5 Variable Binding density</th><th rowspan="2">Unit</th></tr><tr><th>Fe</th><th>Co</th><th>Ni</th><th>Biomass</th><th>Inert</th><th>Precipitate</th></tr></thead><tbody><tr><td>1</td><td>S<sub>su</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>2</td><td>S<sub>aa</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>3</td><td>S<sub>ta</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>4</td><td>S<sub>va</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>5</td><td>S<sub>bu</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>6</td><td>S<sub>pro</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>7</td><td>S<sub>ac</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>8</td><td>S<sub>h2</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>9</td><td>S<sub>ch4</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>10</td><td>S<sub>co2</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>11</td><td>S<sub>h3</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>12</td><td>S<sub>1</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>13</td><td>X<sub>c</sub></td><td>Variable<sup>a</sup></td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>g COD L<sup>-1</sup></td></tr><tr><td>14</td><td>X<sub>ch</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>15</td><td>X<sub>pr</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>16</td><td>X<sub>h1</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>17</td><td>X<sub>su</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>18</td><td>X<sub>aa</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>19</td><td>X<sub>ta</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>20</td><td>X<sub>ca</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>21</td><td>X<sub>pro</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>22</td><td>X<sub>ac</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>23</td><td>X<sub>h2</sub></td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>0.12</td><td>g COD L<sup>-1</sup></td></tr><tr><td>24</td><td>X<sub>1</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>31</td><td>S<sub>hva</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>32</td><td>S<sub>hbu</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>33</td><td>S<sub>hpro</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>34</td><td>S<sub>hac</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>g COD L<sup>-1</sup></td></tr><tr><td>35</td><td>S<sub>hco3</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>M</td></tr><tr><td>36</td><td>S<sub>nh4</sub></td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>M</td></tr><tr><td>37</td><td>S<sub>co32</sub></td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>0.004</td><td>M</td></tr><tr><td>38</td><td>S<sub>Ca2</sub></td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>Variable<sup>a</sup></td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>M</td></tr><tr><td>40</td><td>S<sub>Mg2</sub></td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>Variable<sup>a</sup></td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>0.0025</td><td>M</td></tr><tr><td>42</td><td>S<sub>Ni2</sub></td><td>0.00001</td><td>0.00001</td><td>0.00001</td><td>Variable<sup>a</sup></td><td>0.00001</td><td>0.00001</td><td>0.00001</td><td>0.00001</td><td>0.00001</td><td>M</td></tr><tr><td>44</td><td>S<sub>Co2</sub></td><td>0.0000001</td><td>0.0000001</td><td>Variable<sup>a</sup></td><td>0.0000001</td><td>0.0000001</td><td>0.0000001</td><td>0.0000001</td><td>0.0000001</td><td>0.0000001</td><td>M</td></tr><tr><td>46</td><td>S<sub>Fe2</sub></td><td>0.0001</td><td>Variable<sup>a</sup></td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>M</td></tr><tr><td>51</td><td>S<sub>Fe43</sub></td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>0.006</td><td>M</td></tr><tr><td>58</td><td>S<sub>Is</sub></td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>M</td></tr><tr><td>71</td><td>S<sub>edta4</sub></td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>Variable<sup>a</sup></td><td>0.0001</td><td>0.0001</td><td>0.0001</td><td>M</td></tr></tbody></table>

**Gas transfer processes affecting TEs geochemistry.** The geochemistry of TEs in liquid phase is strongly affected by the H₂S liquid-gas transfer process. The H₂S stripping to gas phase is regulated by the equilibrium concentration of sulfur compounds in the liquid phase. TE interaction with sulfur system during precipitation and adsorption processes drags the equilibrium in the opposite direction towards solid phase. The model also considers liquid-gas transfer processes (Eq. 3) for gaseous components formed during the biological degradation of organic complex matter, namely H₂, CH₄ and CO₂.

**TE speciation during AD.** Modelling the dynamic behavior of TEs during AD is challenging but is necessary to effectively represent the geochemistry of TEs in anaerobic environments. The TE interaction with other compounds (sulfur, phosphorus, carbon, chelating agents and available solid surfaces), and bio-uptake require a substantial investment in model development and verification at different stages. Likewise, TEs interactions with alkali earth elements (Ca/Mg) and other compounds need to be considered, which increase model and computational complexity. Therefore, the speciation model has been developed and tested in stages.
<table><caption>**Table 3.** Parameters used in model comparison with real case scenario.</caption><thead><tr><td>Parameter</td><td>Unit</td><td>Value</td><td>Source</td></tr></thead><tbody><tr><td>Total digester volume</td><td>L</td><td>4</td><td>Literature<sup>43</sup></td></tr><tr><td>Headspace volume</td><td>L</td><td>1</td><td>Literature<sup>43</sup></td></tr><tr><td>OLR</td><td>g VSl<sup>-1</sup> day<sup>-1</sup></td><td>1.45</td><td>Literature<sup>43</sup></td></tr><tr><td>TS</td><td>%</td><td>28.1</td><td>Literature<sup>43</sup></td></tr><tr><td>VS</td><td>% of TS</td><td>95.5</td><td>Literature<sup>43</sup></td></tr><tr><td>TKN</td><td>% of TS</td><td>3.77</td><td>Literature<sup>43</sup></td></tr><tr><td>Total COD</td><td>mg g<sup>-1</sup></td><td>450</td><td>Literature<sup>43</sup></td></tr><tr><td>Temperature</td><td>°C</td><td>37</td><td>Literature<sup>43</sup></td></tr><tr><td>S content of FW</td><td>mol gCOD<sup>-1</sup></td><td>0.0000006</td><td>In this study</td></tr><tr><td>P content of FW</td><td>mol gCOD<sup>-1</sup></td><td>0.00000006</td><td>In this study</td></tr><tr><td>TE (Ni, Co, and Fe) dosed</td><td>mol gCOD<sup>-1</sup></td><td>0 (starvation)</td><td>In this Study</td></tr></tbody></table>

<table><caption>**Table 4.** Various adsorption/desorption reactions considered in the model. *X<sub>biomass</sub> represents the different particulate biomass species considered in the original ADM1.</caption><thead><tr><td>#</td><td>Description</td><td>Reaction</td></tr></thead><tbody><tr><td>1</td><td>Ni-biomass adsorption/desorption</td><td>Ni<sup>2+</sup> + X<sup>*</sup><sub>biomass</sub> ↔ Ni ≡ X<sub>biomass</sub></td></tr><tr><td>2</td><td>Co-biomass adsorption/desorption</td><td>Co<sup>2+</sup> + X<sub>biomass</sub> ↔ Co ≡ X<sub>biomass</sub></td></tr><tr><td>3</td><td>Fe-Biomass Adsorption/Desorption</td><td>Fe<sup>2+</sup> + X<sub>biomass</sub> ↔ Fe ≡ X<sub>biomass</sub></td></tr><tr><td>4</td><td>Ca-biomass adsorption/desorption</td><td>Ca<sup>2+</sup> + X<sub>biomass</sub> ↔ Fe ≡ X<sub>biomass</sub></td></tr><tr><td>5</td><td>Mg-biomass adsorption/desorption</td><td>Mg<sup>2+</sup> + X<sub>biomass</sub> ↔ Fe ≡ X<sub>biomass</sub></td></tr><tr><td>6</td><td>Ni-inert adsorption/desorption</td><td>Ni<sup>2+</sup> + X<sub>inert</sub> ↔ Ni ≡ X<sub>inert</sub></td></tr><tr><td>7</td><td>Co-inert adsorption/desorption</td><td>Co<sup>2+</sup> + X<sub>inert</sub> ↔ Co ≡ X<sub>inert</sub></td></tr><tr><td>8</td><td>Fe-inert adsorption/desorption</td><td>Fe<sup>2+</sup> + X<sub>inert</sub> ↔ Ca ≡ X<sub>inert</sub></td></tr><tr><td>9</td><td>Ca-inert adsorption/desorption</td><td>Ca<sup>2+</sup> + X<sub>inert</sub> ↔ Ca ≡ X<sub>inert</sub></td></tr><tr><td>10</td><td>Mg-inert adsorption/desorption</td><td>Mg<sup>2+</sup> + X<sub>inert</sub> ↔ Mg ≡ X<sub>inert</sub></td></tr><tr><td>11</td><td>Ni-FeS adsorption/desorption</td><td>Ni<sup>2+</sup> + X<sub>FeS</sub> ↔ Ni ≡ X<sub>FeS</sub></td></tr><tr><td>12</td><td>Co-FeS adsorption/desorption</td><td>Co<sup>2+</sup> + X<sub>FeS</sub> ↔ Co ≡ X<sub>FeS</sub></td></tr></tbody></table>

**TE precipitation/dissolution model.** The precipitation model<sup>26</sup> integrates the ADM1 with liquid-solid precipitation reactions. To this aim, new inorganic components and new precipitation/dissolution processes have been incorporated in the ADM1 framework. These new inorganic components influence the biochemical processes by affecting the pH of the system. Although the total free TEs concentration is not directly related to pH, precipitation affects the total free TEs present in an anaerobic digester. New chemical equilibrium association/dissociation, ion pairing, and precipitation reactions have been added to the ADM1. The components of the three chemical systems (carbonate, phosphate and sulfide) react, in the liquid phase, to form precipitates in solid phase. Precipitation has been modeled kinetically by considering a second order law based on the concentration of the participating components and reflecting the important steps of: development of supersaturation, nucleation and growth. Because precipitation reactions are reversible in nature, the dissolution of formed precipitate has been considered. Such dissolution reaction takes place when the system is undersaturated for a particular species and it takes place until the system remains in this condition. The precipitation/dissolution reactions are, thus, governed by the solubility product (K<sub>sp</sub>) values which have been collected from literature and geochemical databases. The full expression of the precipitation/dissolution reaction rates can be found elsewhere<sup>26</sup>.

**TE complexation model.** The complexation module focuses on TE aqueous complexation reactions taking place in AD<sup>37</sup>. In this model, ADM1 has been further modified in order to simulate the TE complexation, precipitation and its effect on AD. In this regard, the biodegradation processes, including hydrolysis, acidogenesis, acetogenesis, methanogenesis, and the acid-base equilibrium reactions have been coupled to the TE complexation, precipitation and redissolution processes. The incorporation of the precipitation and complexation reactions in the physicochemical module led to the definition of new inorganic components (in addition to the precipitation model) in the ADM1 framework. In particular, the complexation of TEs with EDTA and VFAs has been considered. To this aim, EDTA protonation and deprotonation reactions as well as complexation with TEs have been implemented through first order laws. Further, complexation of TEs with VFAs has been considered in the model. The TE-VFA complexes considered in the model are the most likely to be formed in a similar multispecies geochemical environment. The complexation model follows an established and tested modelling approach<sup>35</sup>.
For instance, the overall formation constant of a particular TE-EDTA species is set at a particular value. Consequently, the reverse rate constant is calculated from the stability constant using the following relation:

$$K_{[MeEDTA]}^{2-} = \frac{k_1}{k_{-1}} \qquad (5)$$

Subsequently, the charge balance has been modified in order to consider the effects of new components on the pH of the system. The full expression of the complexation reaction rates can be found elsewhere³⁷.

**TE sorption/desorption model.** TE sorption/desorption has been modeled using a kinetic approach. Specifically, a reversible second order kinetic model has been considered for sorption (forward) and desorption (backward) reactions. Three types of adsorption surfaces have been considered: (1) Biomass; (2) Inerts and (3) Precipitates. Concerning the mineral precipitates, FeS has been recognized as the component playing the most crucial role in sorption processes⁴⁷,⁴⁸. Inert is a major source of sorbent because it harbours many complex polymeric substances, such as lignocellulose, cellulose and various complex protein moieties. The inert is composed of both organic and inorganic material. However, in the proposed model there is no distinction made on this basis. Biomass surfaces can vary widely (carboxyl, hydroxyl, sulfhydryl and phosphoryl functional groups) among different species. To rather complicate, inter species microbial diversity may be recalled. In the proposed model, universal biomass adsorption behaviour has been adopted. It has been hypothesized that the TE bacterial surface interactions are common over all the trophic groups considered in the original ADM1. This has been supported by many studies⁴⁹,⁵⁰.

Seven biomass species, as originally formulated in the ADM1, are construed to participate in the sorption/desorption reactions with all the five metals considered (TEs: Fe, Ni and Co; alkali earth metals: Ca and Mg, Table 4). Inert have been construed to take part in sorption/desorption reactions with Fe, Ni, Co, Mg and Ca, while, only Co and Ni have been considered to participate with precipitates (FeS)⁴⁷,⁴⁸. FeS has been considered as the only mineral precipitate on which adsorption of other TEs takes place which is due to unavailability of supportive studies to incorporate other minerals as adsorbent. Thus, arriving at a system which consists of 35 sorption/desorption reactions for biomass, 5 sorption/desorption reactions for inerts and 2 for sorption/desorption reactions on FeS. The rationale here is to design and implement a mechanistic model which relies on and further extends the ADM1 framework for identifying and selecting important soluble and particulate components for describing TE dynamics.

**Binding sites.** The adsorption module has been conceptualized on: (1) formation of free binding sites; (2) reaction of free binding sites with TEs and alkali earth metals to form occupied binding sites-sorption reaction; (3) desorption of occupied binding sites; (4) competition among TEs for free binding sites (Fig. 8). In the sorption/desorption reactions, a pivoting role is played by the concentration of binding sites⁵¹,⁵². The latter can be found in two states: free, when the sorption reaction has not occurred yet, and occupied when the bound between the adsorbent and adsorbing components is established. The free binding site concentration is found to increase due to the formation of new sorption surfaces such as biomass, inerts and mineral precipitates. Available free TE ions react with these free binding sites leading to the formation of occupied binding sites. The extent of this conversion depends on the stability constant for the particular adsorption species considered. The adsorption process is reversible in nature. The desorption reaction leads to the release of the adsorbed TEs which decay to the bulk pool of free TEs and the formation of free binding sites. The biomass decay leads to the formation of inerts with a concomitant release of TEs adsorbed on the biomass surface. It should be noted that such decay processes are different from desorption processes. Each sorbent component is characterized by a “binding site density”. The binding site density may be defined as the moles of binding sites per unit gram of the respective surface expressed as moles/gram. Unless otherwise mentioned, a binding site density of 2 × 10⁻³² has been used as a default value for biomass and inert surfaces whereas a value of 2 × 10⁻⁵ has been set for precipitate surfaces. It is based on the rational that the density of organic binding sites is higher as compared to inorganic binding sites. The competition of TEs and alkali earth metals for free binding sites has been included in the module in form of competing mass fluxes for particular free binding sites.

**Sorption/desorption kinetics.** The sorption/desorption processes are kinetically controlled and can be represented as follows:

$$TE^{2+} + TES \rightarrow k_1 (TE \equiv TES), \qquad (6)$$

$$(TE \equiv TES) \rightarrow k_{-1} TE^{2+} + TES, \qquad (7)$$

where, k₁ is the TE ≡ TES formation rate constant; k⁻₁ is the dissociation constant. Two pairs of forward and backward reactions in (6 and 7) have equilibrium reactions and constants defined by:

$$TE^{2+} + TES \underset{\stackrel{K_{TE-TES}}{\rightleftharpoons}}^{k_{TE-TE}} (TE \equiv TES), \qquad (8)$$

where, K<sub>TE-TE</sub> is the equilibrium constant for TE-TES. The rate equation for the mechanism in Eq. (8) and introduced in Eq. (1) ρ<sub>s</sub> = ρ<sub>ads</sub> - ρ<sub>des</sub> can be written as:

$$\rho_{ads} = K_{a/d} \cdot k_{ads} \cdot X_{ads} \cdot S_{TE^{2+}}, \qquad (9a)$$
$$ \rho_{des} = K_{a/d} \cdot k_{des} \cdot \bar{X}_{ads}, \qquad (9b) $$

where, $S_{TE^{2+}}$, $X_{ads}$ and $\bar{X}_{ads}$ are the dynamic state variable for free TE concentration, concentration of free binding sites, concentration of occupied binding sites, respectively. $k_{ads}$ is the sorption and $k_{des}$ is the desorption kinetic rate constant. Note that $X_{ads}$ and $\bar{X}_{ads}$ constitute additional state variables of the system, whose dynamics are explicitly tracked through the sorption/desorption rate. The current model formulation considers the biomass, inerts and precipitates as having similar binding characteristics. The model does not distinguish between binding sites based on chemical characteristics. The idea here is to adapt a universal behaviour for inerts, biomass and precipitates based on binding site density which can be latter improved to incorporate specific binding characteristics, e.g. for carboxylate, phosphate and sulfhydryl surfaces.

**Effect of TEs on biochemical processes.** Since TEs have been reported as micronutrients in AD systems and are constituents of co-factors in enzyme systems, their effect on biochemical processes has been explicitly taken into account by introducing an additional term $I_{TE^{2+}}$ in the inhibition expressions used in ADM1$^{53}$. The function considers both the stimulating and inhibiting action of TEs on biochemical rates depending on TE concentrations. This includes an additional non-competitive biostatic inhibition function $I_{TE^{2+}}$, which keeps in track the effect of TEs on the biochemical rates and acts as a growth limiting factor. $I_{TE^{2+}}$ is expressed as:

$$ I_{TE^{2+}} = \frac{a_1 \cdot (TE^{2+} + C^{TE^{2+}}) + a_2}{(TE^{2+} + C^{TE^{2+}})^2 + b_1 \cdot (TE^{2+} + C^{TE^{2+}}) + b_2}, \qquad (10) $$

where $a_1, a_2, b_1, b_2$ are assumed constants which can be adjusted to obtain a desirable optimum dose-response function at a particular TE concentration, $TE^{2+}$ denotes the concentration of bioavailable TEs within the bulk liquid and $C^{TE^{2+}}$ is the fraction of TE which undergoes complexation reaction to form [TE-EDTA]$^{2,54,55}$. Note that the amount of TEs adsorbed on the various surfaces present in AD is not involved in the function $I_{TE^{2+}}$ as the adsorbed TEs are supposed to be unavailable for the biochemical reactions.

**Model parameters and initial conditions.** An in-silico anaerobic batch digester of working volume 0.75 l (head space volume 0.25 l) was selected for the numerical investigations in “Scenario 1: Effect of adsorbents on TE dynamics in anaerobic digestion” to “Scenario 5: Effect of binding site density on TE dynamics in anaerobic digestion”. The temperature of the reactor was set at 35 °C. The initial conditions for the particular scenario were selected on the basis of experience and literature$^{9,22,53}$ which have been reported in Tables 1, 2. The TEs concentration range was selected from literature$^{9}$. The initial values for VFAs were set to zero. The initial values for ammonia and bicarbonate were adjusted at certain value to run the reactor around neutral pH. The initial amount of calcium and magnesium were taken from literature$^{22}$. The values for the saturation constants for sorption processes were assumed in this study. The current model version assumes a single saturation constant for all the sorption/desorption reaction. The values of model parameters for biodegradation, precipitation/dissolution and complexation reactions were taken from literature$^{26,35,37,42,53,56}$.

For model comparison with a real case data$^{43}$, a continuous AD system which studied the effect of TEs on methane production was chosen. The initial conditions for the model simulations were selected by running the reactor for some days with the same OLR. The influent characteristics and operational parameters were obtained from literature$^{43}$ and reported in Table 3. It was assumed that the type and number of reactions for TEs considered in the model remained unaffected with change in the mode of operation. The saturation constants for the adsorption/desorption reactions were kept at the same values as for batch simulations. The binding site density for biomass, release of TE from organic matter and uptake on the biomass were adjusted to suit the experimental behavior of methane production. All the other parameters were kept at the same values as for the batch simulations. Differently from the batch numerical investigations, the inlet and initial EDTA concentrations have been set to zero as the reactor was not augmented with any chelating agent.

**Model implementation.** The system of ordinary differential equations constituting the model has been implemented in an original code and has been solved using the algorithm ODE 15 s, a multistep, variable-order solver on MATLAB platform. Numerical simulations of specific scenarios have been performed to test model reliability.

Received: 14 September 2020; Accepted: 19 February 2021

Published online: 05 April 2021

**References**

1. Appels, L. et al. Anaerobic digestion in global bio-energy production: Potential and research challenges. *Renew. Sustain. Energy Rev.* **15**, 4295–4301 (2011).

2. Appels, L., Baeyens, J., Degrève, J. & Dewil, R. Principles and potential of the anaerobic digestion of waste-activated sludge. *Prog. Energy Combust. Sci.* **34**, 755–781 (2008).

3. Nguyen, D., Gadhamshetty, V., Nitayavardhana, S. & Khanal, S. K. Automatic process control in anaerobic digestion technology: A critical review. *Bioresour. Technol.* **193**, 513–522 (2015).

4. Li, Y., Park, S. Y. & Zhu, J. Solid-state anaerobic digestion for methane production from organic waste. *Renew. Sustain. Energy Rev.* **15**, 821–826 (2011).

5. Fermoso, F. G., Bartacek, J., Jansen, S. & Lens, P. N. L. Metal supplementation to UASB bioreactors: From cell-metal interactions to full-scale application. *Sci. Total Environ.* **407**, 3652–3667 (2009).
6. Mancini, G., Papirio, S., Riccardelli, G., Lens, P. N. L. & Esposito, G. Trace elements dosing and alkaline pretreatment in the anaerobic digestion of rice straw. *Bioresour. Technol.* **247**, 897–903 (2018).

7. Strik, D. P. B. T. B., Domnanovich, A. M., Zani, L., Braun, R. & Holubar, P. Prediction of trace compounds in biogas from anaerobic digestion using the MATLAB Neural Network Toolbox. *Environ. Model. Softw.* **20**, 803–810 (2005).

8. Zandvoort, B. M. H., Van Hullebusch, E. D., Fermoso, F. G. & Lens, P. N. L. Trace metals in anaerobic granular sludge reactors: Bioavailability and dosing strategies. *Eng. Life Sci.* **6**, 293–301 (2006).

9. Thanh, P., Ketheesan, B., Zhou, Y. & Stuckey, D. C. Trace metal speciation and bioavailability in anaerobic digestion: A review. *Biotechnol. Adv.* **34**, 122–136 (2015).

10. Choong, Y. Y., Ismail, N., Abdullah, A. Z. & Yhaya, M. F. Impacts of trace element supplementation on the performance of anaerobic digestion process: A critical review. *Bioresour. Technol.* **209**, 369–379 (2016).

11. Papirio, S. Coupling acid pretreatment and dosing of Ni and Se enhances the biomethane potential of hazelnut skin. *J. Clean. Prod.* **262**, 121407 (2020).

12. Fermoso, F. G. et al. Fate of trace metals in anaerobic digestion. In *Biogas Science and Technology, Advances in Biochemical Engineering/Biotechnology* (eds. Gubitz, G. et al.) **171–195** (2015).

13. Mudhoo, A. & Kumar, S. Effects of heavy metals as stress factors on anaerobic digestion processes and biogas production from biomass. *Int. J. Environ. Sci. Technol.* **10**, 1383–1398 (2013).

14. Ferry, J. G. The chemical biology of methanogenesis. *Planet. Space Sci.* **58**, 1775–1783 (2010).

15. Glass, J. B. & Orphan, V. J. Trace metal requirements for microbial enzymes involved in the production and consumption of methane and nitrous oxide. *Front. Microbiol.* **3**, 1–20 (2012).

16. Kumar, A. G., Kirubagaran, R., Nagesh, N., Rubio, J. P. & Sekaran, G. Influence of fermentation metabolites on redox potential in anaerobic digestion of proteinaceous solid wastes by *Synergistes* sp. *Eng. Life Sci.* **10**, 368–373 (2010).

17. Finazzo, C. et al. Coenzyme B induced coordination of coenzyme M via its thiol group to Ni(I) of F430 in active methyl-coenzyme M reductase. *J. Am. Chem. Soc.* **125**, 4988–4989 (2003).

18. Tallant, T. C. & Krzycki, J. A. Methylthiol:coenzyme M methyltransferase from *Methanosarcina barkeri*, an enzyme of methanogenesis from dimethylsulfide and methylmercaptopropionate. *J. Bacteriol.* **179**, 6902–6911 (1997).

19. Worms, I., Simon, D. F., Hassler, C. S. & Wilkinson, K. J. Bioavailability of trace metals to aquatic microorganisms: importance of chemical, biological and physical processes on biouptake. *Biochimie* **88**, 1721–1731 (2006).

20. Yekta, S. S., Skullberg, U., Danielsson, Å., Björn, A. & Svensson, B. H. Chemical speciation of sulfur and metals in biogas reactors—Implications for cobalt and nickel bio-uptake processes. *J. Hazard. Mater.* **324**, 110–116 (2016).

21. van Hullebusch, E. D. et al. Methodological approaches for fractionation and speciation to estimate trace element bioavailability in engineered anaerobic digestion ecosystems: An overview. *Crit. Rev. Environ. Sci. Technol.* **46**, 1324–1366 (2016).

22. Musvoto, E. V., Wentzel, M. C. M. & Ekama, G. A. M. Integrated chemical-physical processes modelling-II. Simulating aeration treatment of anaerobic digester supernatants. *Water Res.* **34**, 1868–1880 (2000).

23. Batstone, D. J. & Keller, J. Industrial applications of the IWA anaerobic digestion model No.1 (ADM1). *Water Sci. Technol.* **47**, 199–206 (2003).

24. van Rensburg, P., Musvoto, E. V., Wentzel, M. C. & Ekama, G. A. Modelling multiple mineral precipitation in anaerobic digester liquor. *Water Res.* **37**, 3087–3097 (2003).

25. Flores-Alsina, X. et al. Modelling phosphorus (P), sulfur (S) and iron (Fe) interactions for dynamic simulations of anaerobic digestion processes. *Water Res.* **95**, 370–382 (2016).

26. Maharaj, B. C., Mattei, M. R., Frunzo, L., van Hullebusch, E. D. & Esposito, G. ADM1 based mathematical model of trace element precipitation/dissolution in anaerobic digestion processes. *Bioresour. Technol.* **267**, 666–676 (2018).

27. Zhang, Y., Piccard, S. & Zhou, W. Improved ADM1 model for anaerobic digestion process considering physico-chemical reactions. *Bioresour. Technol.* **196**, 279–289 (2015).

28. Huber, P., Neyret, C. & Fourest, E. Implementation of the anaerobic digestion model (ADM1) in the PHREEQC chemistry engine. *Water Sci. Technol.* **76**, 1090–1103 (2017).

29. Maharaj, B. C. et al. Mathematical modelling of trace element dynamics in anaerobic digestion environments. In *Trace Elements in Anaerobic Biotechnologies* 101–152 (IWA Publishing, 2019).

30. Vaneeckhaute, C. et al. Development, implementation, and validation of a generic nutrient recovery model (NRM) library. *Environ. Model. Softw.* **99**, 170–209 (2018).

31. Laera, A. et al. A simultaneous assessment of organic matter and trace elements bio-accessibility in substrate and digestate from an anaerobic digestion plant. *Bioresour. Technol.* **288**, 121587 (2019).

32. Ortner, M., Rachbauer, L., Somitsch, W. & Fuchs, W. Can bioavailability of trace nutrients be measured in anaerobic digestion?. *Appl. Energy* **126**, 190–198 (2014).

33. Batstone, D. J. et al. Towards a generalized physicochemical framework. *Water Sci. Technol.* **66**, 1147–1161 (2012).

34. Frunzo, L., Fermoso, F. G., Luongo, V., Mattei, M. R. & Esposito, G. ADM1-based mechanistic model for the role of trace elements in anaerobic digestion processes. *J. Environ. Manag.* **241**, 587–602 (2019).

35. Willet, A. I. & Rittmann, B. E. Slow complexation kinetics for ferric iron and EDTA complexes make EDTA non-biodegradable. *Biodegradation* **14**, 105–121 (2003).

36. Hauduc, H., Smith, S., Szabo, A., Murthy, S. & Tak, I. A dynamic physicochemical model for chemical phosphorus removal. *Water Res.* **73**, 157–170 (2015).

37. Maharaj, B. C., Mattei, M. R., Frunzo, L., Hullebusch, E. D. V. & Esposito, G. ADM1 based mathematical model of trace element complexation in anaerobic digestion processes. *Bioresour. Technol.* **276**, 253–259 (2019).

38. Kazadi Mbamba, C., Batstone, D. J., Flores-Alsina, X. & Tait, S. A generalised chemical precipitation modelling approach in wastewater treatment applied to calcite. *Water Res.* **68**, 342–353 (2015).

39. Bartacek, J., Fermoso, F. G., Catena, A. B. & Lens, P. N. L. Effect of sorption kinetics on nickel toxicity in methanogenic granular sludge. *J. Hazard. Mater.* **180**, 289–296 (2010).

40. Jong, T. & Parry, D. L. Adsorption of Pb (II), Cu (II), Cd (II), Zn (II), Ni (II), Fe (II), and As (V) on bacterially produced metal sulfides. *J. Colloid Interface Sci.* **275**, 61–71 (2004).

41. Thanh, P., Ketheesan, B., Zhou, Y. & Stuckey, D. C. Effect of ethylenediamine-N,N'-disuccinic acid (EDDS) on the speciation and bioavailability of Fe2+ in the presence of sulfide in anaerobic digestion. *Bioresour. Technol.* **229**, 169–179 (2017).

42. Schwarz, A. O. & Rittmann, B. E. A biogeochemical framework for metal detoxification in sulfidic systems. *Biodegradation* **18**, 675–692 (2007).

43. Climenhaga, M. A. & Banks, C. J. Anaerobic digestion of catering wastes: Effect of micronutrients and retention time. *Water Sci. Technol.* **57**, 687–693 (2008).

44. Grass, G., Fricke, B. & Nies, D. H. Control of expression of a periplasmic nickel efflux pump by periplasmic nickel concentrations. *Biometals* **18**, 437–448 (2005).

45. Zhang, W., Zhang, L. & Li, A. Enhanced anaerobic digestion of food waste by trace metal elements supplementation and reduced metals dosage by green chelating agent [S,S]-EDDS via improving metals bioavailability. *Water Res.* **84**, 266–277 (2015).

46. Race, M. Applicability of alkaline precipitation for the recovery of EDDS spent solution. *J. Environ. Manag.* **203**, 358–363 (2017).

47. van Hullebusch, E. D. et al. Cobalt sorption onto anaerobic granular sludge: Isotherm and spatial localization analysis. *J. Biotechnol.* **121**, 227–240 (2006).
48. Watson, J. H. P. et al. Heavy metal adsorption on bacterially produced FeS. *Miner. Eng.* **8**, 1097–1108 (1995).

49. Jiang, B. *The Effect of Trace Elements on the Metabolism of Methanogenic Consortia* (Wageningen University, 2006).

50. Johnson, K. J., Szymanowski, J. E. S., Borrok, D., Huynh, T. Q. & Fein, J. B. Proton and metal adsorption onto bacterial consortia: Stability constants for metal—bacterial surface complexes. *Chem. Geol.* **239**, 13–26 (2007).

51. D'Acunto, B., Frunzo, L., Luongo, V. & Mattei, M. R. Modeling heavy metal sorption and interaction in a multispecies biofilm. *Mathematics* **7**, 781 (2019).

52. D'Acunto, B., Frunzo, L. & Mattei, M. R. On a free boundary problem for biosorption in biofilms. *Nonlinear Anal. Real World Appl.* **39**, 120–141 (2018).

53. Batstone, D. J. et al. The IWA anaerobic digestion model no 1 (ADM1). *Water Sci. Technol.* **45**, 65–73 (2002).

54. Vintiloiu, A. et al. Effect of ethylenediaminetetraacetic acid (EDTA) on the bioavailability of trace elements during anaerobic digestion. *Chem. Eng. J.* **223**, 436–441 (2013).

55. Hu, Q. H., Li, X. F., Liu, H., Du, G. C. & Chen, J. Enhancement of methane fermentation in the presence of Ni²⁺ chelators. *Biochem. Eng. J.* **38**, 98–104 (2008).

56. Callander, I. J. & Barford, J. P. Precipitation, chelation, and the availability of metals as nutrients in anaerobic digestion. II. Applications. *Biotechnol. Bioeng.* XXV, 1959–1972 (1983).

## Acknowledgements

This study was supported by the Marie Skłodowska-Curie European Joint Doctorate-ABWET, funded by the Horizon 2020 program [Grant agreement 643071]. The authors thank: *COST Action* 1302 (European Network on Ecological Roles of Trace *Metals* in Anaerobic Biotechnologies), Dr. Ludovico Pontoni and Dr. Vincenzo Luongo for engaging discussions.

## Author contributions

B.C.M. carried out the mathematical model formulation, identification of the variables and processes involved, the definition of the modelling scenarios, implementation of the numerical code, verification (chemical equilibrium calculations), interpretation of simulations results, and wrote the manuscript. M.R.M. and L.F. carried out conceptualization and designing of the model framework, implementation of the numerical code, interpretation of the simulation results and thoroughly revised the manuscript. G.E. and E.vH. participated in the identification of the variables and processes involved, in the definition of the modelling scenarios, in the results interpretation and thoroughly revised the manuscript. G.E. organized and coordinated the research funding.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to B.C.M.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021