A Dynamic Biofilm Model for a Microbial Electrolysis Cell

Abstract

In this work, a mathematical description of a Microbial Electrolysis Cell (MEC) is proposed, taking into account the global mass balances of the different species in the system and considering that all the involved microorganisms are attached to the anodic biological film. Three main biological reactions are introduced, which were obtained from the solution of partial differential equations describing the spatial distribution of potential and substrate in the biofilm. The simulation of the model was carried out using numerical methods, and the results are discussed.

1. Introduction

Bioelectrochemical Systems (BES) technology has been of significant interest recently due to the advantage of producing energy from renewable organic materials. Mainly, Microbial Electrochemical Cells (MXC) have gained attention among BES due to its capacity to treat wastewater and simultaneously recover energy [1,2].

MXC processes take advantage of Anode-Respiring Bacteria (ARB) applied in electrochemical systems. During the respiration process, an electrical current is produced, and it is used in one way or another for producing energy depending on the type of MXC system. Two processes among them have received significant attention from many researchers in the last two decades due to the technological and economical advantages obtained from electricity generation in the case of Microbial Fuel Cells (MFC) and hydrogen production in the case of Microbial Electrolysis Cells (MEC). The current produced in the MFC system is used for the reduction of oxidized compounds (mainly oxygen); in MEC systems, the current is used to produce hydrogen gas [3,4]. In this sense, the MEC process has been recently considered as a promising renewable energy alternative to fossil fuels.

An MEC is an electrochemical system that generates hydrogen and removes organic matter from wastewater with the help of an externally-applied potential and the activity of exoelectrogens microorganisms. The microorganisms oxidize the organic matter and transfer the electrons outside their cells as a part of their respiration process. The released electrons go to the anode, where the applied potential increases their energy level. Simultaneously, protons are released at the anode and transferred to the cathode through the bulk phase or a membrane (if any). Finally, the electrons carry out the reduction of protons at the cathode. Due to the activity of microorganisms, in theory, it is necessary to apply externally only 135 mV to produce the hydrogen, which is considerably lower than the 1210 mV thermodynamically required by the classical electrolysis of water [5,6]. Hence, the advantages of MEC are the lower voltage needed for the H₂ generation and the removal of organic matter from wastewater.

The performance of BES depends on several electrochemical, physical-chemical, and biological factors. The understanding of the dynamic relations among these is critical to making this technology more efficient [7,8,9]. Furthermore, the mathematical description of these relations is fundamental to the design and optimization of this kind of process [10,11]. Indeed, modeling of MXC processes has been an active area throughout the last decade. Recently, a broad classification of mathematical models of BES systems has been reported by Gadkari et al. [12]. Nevertheless, most contributions are mainly focused on MFC, while modeling research work on MEC systems is still limited [10,12]. In any case, differential equations are used as the main mathematical tools for describing MXC systems.

Based on a previously-reported MFC model, Pinto et al. [13] described a MEC system through a set of eight Ordinary Differential Equations (ODE) [14]. The mathematical description includes: (i) a microbial multi-population dynamic (anodophilic, fermentative, methanogenic acetoclastic, and methanogenic hydrogenophilic microorganisms); (ii) a simplification of anaerobic degradation of wastewater based on two single steps (hydrolysis and fermentation); (iii) Bernard’s anaerobic digestion kinetics [15]; and (iv) hydrogen consumption and production. Afterward, the MEC model was modified with additional considerations intended for a unified description of both MFC and MEC systems [16]. In both models, there is not a dynamical description of biofilm.

On the other hand, based on a previously-reported MFC model [17], Alavijeh et al. [18] described a one-dimensional MEC system through a set of Partial Differential Equations (PDE). The model considers several variables of liquid bulk and biofilm, such as: (i) a microbial multi-population dynamic (anodophilic, fermentative, and methanogenic acetoclastic microorganisms); (ii) active and inactive microorganisms; (iii) spatial one-dimensional description of substrates and local potential in the biofilm; (iv) a one-dimensional dynamical description of biofilm thickness; (v) global dynamical description of substrates in the liquid bulk. The model provides a detailed description of the MEC system. Nevertheless, despite the model being simplified with linear boundary conditions, the numerical solution involves solving the nonlinear, time-dependent, and coupled system of partial and ordinary differential equations. Then, the numerical solution requires numerical in-out cycles among several numerical methods (trust-region-reflective method, Runge–Kutta–Fehlberg method, and Gear’s method).

After that, Mardanpour and Yaghmaei [19] proposed a dynamic mathematical model to study microfluidic MEC systems with potential application in biohydrogen production for implantable medical devices. The mathematical description includes: (i) one microbial population dynamic; (ii) active and inactive microorganisms; (iii) one-dimensional local electrical potential description. Despite a detailed description of microbial biomass particles, there is not a differential equation related to biofilm thickness; instead, an algebraic equation for biofilm thickness is presented. The model provides an accurate description of the dynamical MEC system; however, it is limited to microfluidic applications.

Recently, Hussain et al. [20] described the MEC system as a simple MEC Equivalent Electrical Circuit (EEC) model. Two first order differential equations were proposed to represent the voltage dynamics of internal capacitors of the EEC. The main advantage of this proposal is the online procedure for estimating EEC model parameters. Indeed, the approach can be used to achieve real-time monitoring of the MEC system. Despite variation in chemical oxygen demand concentration being captured by the proposed monitoring method, non-information of the microorganism and biofilm dynamics was described.

Mathematical models that take into account biofilm growth dynamics of a mixed microbial population, without the drawback of high computational resources in the model simulation, should be advantageous in the MEC system for optimization and control purposes. However, to the best of our knowledge, a mathematical MEC model based on those ideas has not been reported. In this work, we propose a dynamical model for an MEC with an emphasis on biofilm growth dynamics of a mixed population of microbes and biofilm processes such as the substrate diffusion, the potential loss, and the detachment of biomass. To avoid high computational resources in the MEC model simulation, the dynamical behavior of substrate distribution and potential inside the biofilm (PDE) are included in a set of ODE systems using average reaction rates.

The paper is organized as follows: the general MEC process is described in Section 2. Then, the proposed dynamical model is shown in Section 3. The implementation via simulation runs of the model for two study cases, chronoamperometry and voltammetry, is addressed in Section 4. In Section 5, some potential application and further comments are provided. Finally, some concluding remarks are presented in Section 6.

2. The Cell

In this work, an MEC system composed of two chambers, the anodic and the cathodic one, separated by an ionic membrane is considered. Both of them are assumed under anaerobic conditions to avoid the oxygen being the final electron acceptor.

In the anodic chamber, there is an electrode in which the biomass attaches and creates a microbial film (named the bioanode). The biofilm takes the nutrients from the bulk phase in the anodic chamber and produces the electrons that travel externally to the cathode. Bacteria are capable of transferring electrons outside the cell, either directly or by endogenous mediators (which can also be added), named exoelectrogens [3]. These latter include a wealth of genera such as Geobacter, Shewanella, Pseudomonas, and others [21].

The absence of substances that can be final electron acceptors (e.g., sulfates, nitrates, metal cations, etc.) is assumed. Protons are also bio-products, and they travel to the cathodic chamber through the membrane, commonly Nafion; a proton exchange membrane [5,21,22,23].

The produced electrons pass through an external potential, which increases their energy level, and finally go to the cathode. If they have adequate potential, the protons coming from the solution in the anodic chamber are reduced, and the molecular hydrogen is produced.

Regarding the produced electrons, there is an open question related to the transference of electrons from the microbial cell to the anode. Rozendal et al. [22] stated that this is either by direct contact with an electrode surface or aided by excreted redox mediators. The most recent studies have shown that the mechanisms for the electron transfer rely on the type of microorganism [24].

3. Dynamical Model

The model was obtained from a time-dependent mass balance of the different species in the biofilm and the anodic chamber. In the biofilm, a mass balance was also made for the spatial distribution of substrate and a charge balance for the potential spatial distribution. For the applied potential, the anode was considered the working electrode.

3.1. The System

Each of the two chambers contained a volume of fluid $V_a$ (mL) for the anodic one and $V_c$ (mL) for the cathodic one, respectively. The anode chamber was fed continuously with a volumetric flow F (mL d⁻¹). It was considered that the supplied solution was composed of a phosphate buffer (pH = 7.0) containing oligo-elements [21,22] and acetate as the only carbon source with a concentration $S_{in}$ (mmol ${\mathrm{C}}_2{\mathrm{H}}_3{\mathrm{O}}_2^{-}$ mL⁻¹). The anode had an active area $A_a$ (cm²) where the microorganisms attached and created a uniform biofilm of thickness $L_f$ (cm). The cathodic chamber contained only an abiotic buffered solution (pH = 7.0).

The electrodes were externally connected by a potentiometer that kept the potential $E_a$ (V) vs. the reference of the anode constant. It was assumed that the reference electrode and the anode were close enough that the ohmic losses between them were negligible.

3.2. Biomass

As the only carbon source was acetate, just two populations were considered in the system: the methanogenic and the exoelectrogens microorganisms. Both of them competed for the substrate and space in the biofilm.

Due to the absence of molecular hydrogen in the anodic chamber, the methane was produced by the acetoclastic process as follows:

$${\mathrm{CH}}_3{\mathrm{COO}}^{-}+{\mathrm{H}}_2\mathrm{O}\to {\mathrm{HCO}}_3^{-}+{\mathrm{CH}}_4$$

The exoelectrogens microorganisms carried out the following reaction:

$${\mathrm{CH}}_3{\mathrm{COO}}^{-}+4{\mathrm{H}}_2\mathrm{O}\to 2{\mathrm{HCO}}_3^{-}+9{\mathrm{H}}^{+}+8{\mathrm{e}}^{-}$$

The microorganisms used part of the produced electrons for self-growth and their own energy production [25]; this fraction is denoted as $f_s$, and it was assumed that the rest of the fraction, $f_{exp}=1-f_s$, was expelled directly to the anode.

It was considered that all the microorganisms were attached to the biofilm and there was no activity in the liquid phase. The concentration in the biofilm was $x_a$ (mg VS cm⁻²) for the acetoclastic methanogens and $x_e$ (mg VS cm⁻²) for the exoelectrogens microorganisms, where $VS$ (mg) stands for the measurement of the biomass as volatile solids.

3.3. Biological Reaction Pathways

The considered biological reactions were acetoclastic methanogenesis (with a reaction rate $r_a={\mu }_a\left(S\right)x_a$):

$$\mathrm{S}\stackrel{r_a}{\to }{k}_1{x}_a+k_2{\mathrm{CH}}_4+k_3{\mathrm{CO}}_2.$$

electrogenesis (with a reaction rate $r_e={\mu }_e(S,E_a)x_e$:

$$\mathrm{S}\stackrel{r_e}{\to }{k}_4{x}_e+k_5{\mathrm{e}}^{-}+k_6{\mathrm{CO}}_2.$$

and self-oxidation of biomass through endogenous respiration (with a reaction rate $r_{s-o}=r_{res}\left(E_a\right)x_e$:

$$x_e\stackrel{r_{s-o}}{\to }{k}_7{\mathrm{e}}^{-}+k_8{\mathrm{CO}}_2$$

where $k_i$ ($i=1,2,\dots ,8$) are the respective yield coefficients. In the ideal case, when no biomass was produced and all the substrate was oxidized, their values were obtained from carbon and electron balances. For details, see Appendix A.

3.4. Microorganisms’ Kinetics

The rate at which methanogenic microorganisms consumed the substrate follows the Haldane kinetics [15]:

$${\mu }_a\left(S\right)={\mu }_{a,max}\frac{S}{S+K_{S,a}+\frac{S^2}{K_I}}$$

(1)

where ${\mu }_{a,max}$ is the maximum specific growth rate without inhibition (mmol S mg VS⁻¹ d⁻¹), $K_{S,a}$ is the half maximum rate concentration of S without inhibition (mmol S mL⁻¹), and $K_I$ is the inhibition constant associated with S (mmol S mL⁻¹).

The growth of the exoelectrogens microorganisms was affected by the substrate concentration and the potential of the anode (which can be seen as the electron acceptor). The rate is described by the multiplicative Nernst–Monod kinetics [17,26]:

$${\mu }_e(S,E_a)={\mu }_{e,max}\frac{S}{S+K_{S,e}}\frac{1}{1+\exp \{-\frac{\mathfrak{F}}{RT}(E_a-E_{K_A})\}}$$

(2)

where: ${\mu }_{e,max}$ is the maximum rate of utilization by the exoelectrogens microorganisms (mmol S mg VS⁻¹ d⁻¹), S is the acetate concentration (mmol S mL⁻¹), $K_{S,e}$ is the half-maximum-rate-concentration of species S (mmol S mL⁻¹), $E_a$ (V) is the anode potential, and $E_{K_A}$ (V) is the potential at which ${\mu }_e=\frac{1}{2}{\mu }_{e,max}$ when $S\gg K_{S,e}$. Both potentials are expressed vs. the same reference electrode.

Besides the growth, there is a process by which the microorganisms die or stop working. This process is the inactivation of active biomass, and it was considered that it followed the next kinetics [17,27]:

$$r_{in}=b_{in}(x_a+x_e)$$

(3)

where $r_{in}$ $\left(\mathrm{mgVS}{\mathrm{cm}}^{-2}{\mathrm{d}}^{-1}\right)$ is the inactivation rate, and $b_{in}$ (d⁻¹) is the inactivation constant, and it is assumed to be equal for both microbial populations. Notice that $b_{in}$ can also be defined as a function of the medium characteristics (pH, T, the presence of other substances); hence, $b_{in}$ could be a time-dependent parameter for a batch reactor where the qualities of the medium may change with time. However, for this work, it is considered constant.

The endogenous respiration is related to the self-oxidation of exoelectrogens microorganisms [17], and it is described using the Nernst–Monod equation:

$$r_{res}\left(E_a\right)=b_{res}\frac{1}{1+\exp \{-\frac{\mathfrak{F}}{RT}(E_a-E_{K_A})\}}$$

(4)

where $b_{res}$ (d⁻¹) is the endogenous decay coefficient.

3.5. Biofilm

The biofilm was composed of three different species: the two active microorganisms $x_a$ and $x_e$; and the inert biomass $x_i$ (mg VS cm⁻²). The latter is due to the inactive microorganisms, extracellular polymeric substances, etc.

The total volume (cm³) of the biofilm was:

$$V_{bio}=A_a{L}_f$$

and its mass (mg) was:

$$m_{bio}={\rho }_x{V}_{bio}=A_a(x_a+x_e+x_i)$$

where ${\rho }_x$ is the biofilm density (mg VS cm⁻³); it was assumed to be constant and uniform through the film.

From the previous expressions, the concentration of each species in the biofilm can be calculated as follows:

$$x_a=L_f{\rho }_x{\varphi }_a$$

(5)

$$x_e=L_f{\rho }_x{\varphi }_e$$

(6)

$$x_i=L_f{\rho }_x(1-{\varphi }_a-{\varphi }_e)$$

(7)

where ${\varphi }_a$ and ${\varphi }_e$ are the mass fractions of $x_a$ and $x_e$, respectively.

It was assumed that there was no spatial variation of biomass along the biofilm. However, the substrate concentration and the potential can vary inside the film due to the mass and charge transfers [17]. To find these variations, some balances were made for the acetate and electrons. It was assumed that all the transfers were unidirectional (axis z) and perpendicular to the anode area. The anode surface was defined at $z=0$, and the interface between the biofilm and the liquid phase was located at $z=L_f$.

3.5.1. Substrate Diffusion

As the acetate reacted, a concentration gradient appeared along the biofilm thickness. Departing from the variation in the substrate concentration ($S^{bio}$), the following mass balance can be established:

$$\frac{\partial S^{bio}}{\partial t}=\frac{\partial J_z}{\partial z}+{\rho }_x[{\mu }_a\left(S^{bio}\right){\varphi }_a+{\mu }_e(S^{bio},E_a^{bio}){\varphi }_e]$$

where $J_z$ (mmol S cm⁻² d⁻¹) is the substrate flux in the z-direction, i.e.,

$$J_z=-D_e\frac{\partial S^{bio}}{\partial z}$$

with $D_e$ (cm⁻² d⁻¹) the effective diffusion of the substrate in the biofilm. Due to the small thickness of the film, once the concentration of the substrate in the bulk changed, it was assumed that the distribution of $S^{bio}$ inside the biofilm immediately reached its steady state. Then, the previous equation can be expressed as:

$$D_e\frac{{\partial }^2{S}^{bio}}{\partial z^2}-{\rho }_x[{\mu }_a\left(S^{bio}\right){\varphi }_a+{\mu }_e(S^{bio},E_a^{bio}){\varphi }_e]=0$$

(8)

Notice that the functions ${\mu }_a\left(S^{bio}\right)$ and ${\mu }_e(S^{bio},E_a^{bio})$ are indeed Equations (1) and (2) evaluated on $S^{bio}\left(z\right)$ and $E_a^{bio}\left(z\right)$.

Due to the absence of substrate diffusion to the electrode, the first boundary condition was:

$${\left.J_z\right|}_{z=0}=-D_e{\left.\frac{\partial S^{bio}}{\partial z}\right|}_{z=0}=0$$

It was assumed that there were no concentration gradients between the bulk and the biofilm surface. Then, the second boundary condition was:

$$S^{bio}{|}_{z=L_f}=S$$

where S is the substrate concentration in the bulk of the anodic chamber.

3.5.2. Potential Losses

As the electrons flowed through the biofilm, there was a potential loss. The equation describing the potential variation $E_a^{bio}$ through the biofilm is:

$$\frac{\partial E_a^{bio}}{\partial t}=\frac{\partial j_z}{\partial z}+\frac{\mathfrak{F}}{\gamma }\left[k_5{\mu }_e\left(S^{bio},E_a^{bio}\right)+k_7{r}_{res}\left(E_a^{bio}\right)\right]{\rho }_x{\varphi }_e$$

where $j_z\left(\mathrm{mA}{\mathrm{cm}}^{-2}\right)$ is the flux of electrons (current density) in the z-direction, i.e.,

$$j_z=-{\kappa }_{bio}\frac{\partial E_a^{bio}}{\partial z}$$

Here, it was considered that electrons flowing through the biofilm immediately reached their steady-state. Then, the equation describing the potential variation through the biofilm is:

$${\kappa }_{bio}\frac{{\partial }^2{E}_a^{bio}}{\partial z^2}-\frac{\mathfrak{F}}{\gamma }\left[k_5{\mu }_e\left(S^{bio},E_a^{bio}\right)+k_7{r}_{res}\left(E_a^{bio}\right)\right]{\rho }_x{\varphi }_e=0$$

(9)

where $\mathfrak{F}$ is the Faraday constant (96,487 mC mmol $e^{-1}$); $\gamma$ is the conversion time factor (86,400 s d⁻¹); ${\kappa }_{bio}$ (mS cm⁻¹ = mA V⁻¹ cm⁻¹) is the biofilm conductivity that can be seen as the sum of the effects related to the electrons’ transfer mechanisms.

The boundary conditions related to a fixed potential at the anode were:

$$E_a^{bio}{|}_{z=0}=E_a$$

Furthermore, it was assumed that electrons conducted only on the biofilm matrix; then, a second boundary condition was:

$${\left.j_z\right|}_{z=L_f}=-{\kappa }_{bio}{\left.\frac{\partial E_a^{bio}}{\partial z}\right|}_{z=L_f}=0$$

Thus, the substrate and the potential distribution along the biofilm, $S^{bio}\left(z\right)$ and $E_a^{bio}\left(z\right)$, were obtained by solving Equations (8) and (9).

3.5.3. Detachment

It was assumed that the biofilm thickness was finite [13,28]. Then, there was a mechanism by which biomass was detached from the film and released into the bulk liquid. Thus, it was assumed that the rate of detachment $r_{det}$ $\left(\mathrm{mgVS}{\mathrm{cm}}^{-2}{\mathrm{d}}^{-1}\right)$ followed a second-order dependence on the biofilm thickness as follows [29]:

$$r_{det}=b_{det}{\rho }_x{L}_f^2$$

(10)

where $b_{det}$ is the detachment constant (cm⁻¹ d⁻¹), which can be dependent on the shear stress.

3.6. Mass Balances

3.6.1. Average Reaction Rates

Due to the mass and charge transfers, the active microorganisms along the biofilm were exposed to different substrate concentrations and potentials. Therefore, their growth and oxidation rates depended on the position inside the film. The reaction rates affected by the z-position in the biofilm were:

$${\mu }_a\left(S^{bio}\left(z\right)\right),{\mu }_e(S^{bio}\left(z\right),E_a^{bio}\left(z\right)),r_{res}\left(E_a^{bio}\left(z\right)\right)$$

To find the average reaction rates, the previous expressions can be integrated along the biofilm as follows:

$${\overline{\mu }}_a\left(S^{bio}\right)=\frac{1}{L_f}{\int }_0^{L_f}{\mu }_a\left(S^{bio}\left(z\right)\right)dz$$

(11)

$${\overline{\mu }}_e(S^{bio},E_a^{bio})=\frac{1}{L_f}{\int }_0^{L_f}{\mu }_e(S^{bio}\left(z\right),E_a^{bio}\left(z\right))dz$$

(12)

$${\overline{r}}_{res}\left(E_a^{bio}\right)=\frac{1}{L_f}{\int }_0^{L_f}{r}_{res}\left(E_a^{bio}\left(z\right)\right)dz$$

(13)

Notice that it cannot be expected that the average reaction rates are strictly as accurate as the complete local reaction rate profile; however, the proposed approach intended to capture the entire dynamic related to the reaction rate profile.

3.6.2. Global Balance of Species in the MEC

The dynamical behavior of the system was obtained from a mass balance of the different species: (i) substrate; (ii) biofilm; (iii) active acetoclastic biomass; and (iv) active exoelectrogens microorganisms.

The accumulation of the substrate inside the biofilm was neglected due to the small volume of the biofilm compared to the liquid phase. Then, the mass balance for the substrate can be expressed as:

$$\frac{d}{dt}\left(V_{a}S\right)=F[S^{in}-S]-A_a[{\overline{\mu }}_a\left(S^{bio}\right)x_a+{\overline{\mu }}_e(S^{bio},E_a^{bio})x_e]$$

A constant volume of the anodic chamber was assumed; then, the previous equation can be rearranged as:

$$\frac{dS}{dt}=\frac{F}{V_a}[S^{in}-S]-\frac{A_a}{V_a}[{\overline{\mu }}_a\left(S^{bio}\right)L_f{\rho }_x{\varphi }_a+{\overline{\mu }}_e(S^{bio},E_a^{bio})L_f{\rho }_x{\varphi }_e]$$

(14)

The mass balance for the biofilm follows as:

$$\frac{d}{dt}\left({\rho }_x{A}_a{L}_f\right)=A_a[k_1{\overline{\mu }}_a\left(S^{bio}\right)x_a+k_4{\overline{\mu }}_e(S^{bio},E_a^{bio})x_e-b_{det}{\rho }_x{L}_f^2-{\overline{r}}_{res}\left(E_a^{bio}\right)x_e]$$

$$\frac{d{L}_f}{dt}=k_1{\overline{\mu }}_a\left(S^{bio}\right)L_f{\varphi }_a+k_4{\overline{\mu }}_e(S^{bio},E_a^{bio})L_f{\varphi }_e-b_{det}{L}_f^2-{\overline{r}}_{res}\left(E_a^{bio}\right)L_f{\varphi }_e$$

(15)

The mass balance for the active acetoclastic biomass can be expressed as:

$$\frac{d}{dt}\left(A_a{x}_a\right)=A_a{k}_1{\overline{\mu }}_a\left(S^{bio}\right)x_a-A_a{b}_{in}{x}_a-A_a{r}_{det}{\varphi }_a$$

$$\frac{d}{dt}\left(A_a{L}_f{\rho }_x{\varphi }_a\right)=A_a{k}_1{\overline{\mu }}_a\left(S^{bio}\right)L_f{\rho }_x{\varphi }_a-A_a{b}_{in}{L}_f{\rho }_x{\varphi }_a-A_a{b}_{det}{\rho }_x{L}_f^2{\varphi }_a$$

$$\frac{d}{dt}\left(L_f{\varphi }_a\right)=k_1{\overline{\mu }}_a\left(S^{bio}\right)L_f{\varphi }_a-b_{in}{L}_f{\varphi }_a-b_{det}{L}_f^2{\varphi }_a$$

$$L_f\frac{d{\varphi }_a}{dt}+{\varphi }_a\frac{d{L}_f}{dt}=k_1{\overline{\mu }}_a\left(S^{bio}\right)L_f{\varphi }_a-b_{in}{L}_f{\varphi }_a-b_{det}{L}_f^2{\varphi }_a$$

$$\frac{d{\varphi }_a}{dt}=k_1{\overline{\mu }}_a\left(S^{bio}\right){\varphi }_a-b_{in}{\varphi }_a-b_{det}{L}_f{\varphi }_a-\frac{{\varphi }_a}{L_f}\frac{d{L}_f}{dt}$$

From the substitution of Equation (15) in the previous equation, it is possible to obtain:

$$\frac{d{\varphi }_a}{dt}=k_1{\overline{\mu }}_a\left(S^{bio}\right)[{\varphi }_a-{\varphi }_a^2]-b_{in}{\varphi }_a-k_4{\overline{\mu }}_e(S^{bio},E_a^{bio}){\varphi }_e{\varphi }_a+r_{res}{\varphi }_e{\varphi }_a.$$

(16)

The mass balance for the active exoelectrogens microorganisms can be expressed as:

$$\frac{d}{dt}\left(A_a{x}_e\right)=A_a{k}_4{\overline{\mu }}_e(S^{bio},E_a^{bio})x_e-A_a[b_{in}+{\overline{r}}_{res}\left(E_a^{bio}\right)]x_e-A_a{r}_{det}{\varphi }_e$$

$$\frac{d}{dt}\left(A_a{L}_f{\rho }_x{\varphi }_e\right)=A_a{k}_4{\overline{\mu }}_e(S^{bio},E_a^{bio})L_f{\rho }_x{\varphi }_e-A_a[b_{in}+{\overline{r}}_{res}\left(E_a^{bio}\right)]L_f{\rho }_x{\varphi }_e-A_a{b}_{det}{\rho }_x{L}_f^2{\varphi }_e$$

$$\frac{d}{dt}\left(L_f{\varphi }_e\right)=k_4{\overline{\mu }}_e(S^{bio},E_a^{bio})L_f{\varphi }_e-[b_{in}+{\overline{r}}_{res}\left(E_a^{bio}\right)]L_f{\varphi }_e-b_{det}{L}_f^2{\varphi }_e$$

$$L_f\frac{d{\varphi }_e}{dt}+{\varphi }_e\frac{d{L}_f}{dt}=k_4{\overline{\mu }}_e(S^{bio},E_a^{bio})L_f{\varphi }_e-[b_{in}+{\overline{r}}_{res}\left(E_a^{bio}\right)]L_f{\varphi }_e-b_{det}{L}_f^2{\varphi }_e$$

$$\frac{d{\varphi }_e}{dt}=k_4{\overline{\mu }}_e(S^{bio},E_a^{bio}){\varphi }_e-[b_{in}+{\overline{r}}_{res}\left(E_a^{bio}\right)]{\varphi }_e-b_{det}{L}_f{\varphi }_e-\frac{{\varphi }_e}{L_f}\frac{d{L}_f}{dt}$$

From the substitution of Equation (15) in the previous equation, it is possible to obtain:

$$\frac{d{\varphi }_e}{dt}=[k_4{\overline{\mu }}_e(S^{bio},E_a^{bio})-{\overline{r}}_{res}\left(E_a^{bio}\right)][{\varphi }_e-{\varphi }_e^2]-b_{in}{\varphi }_e-k_1{\overline{\mu }}_a\left(S^{bio}\right){\varphi }_a{\varphi }_e$$

(17)

where the states of the MEC model are the substrate concentration in the liquid phase S (g/L); the thickness of biomass in the film $L_f$ (m); and the mass fractions of microorganisms ${\varphi }_a$ and ${\varphi }_e$. The average reaction rates were obtained from the Equations (11)–(13). Then, their calculation required the integration of Equations (8) and (9) for each instant of time. For details about the integration method to solve this set of ordinary differential equations, see Appendix B.

3.7. Current

Under the assumption that all the generated electrons go to the anode, the expected current I (mA) was:

$$I=\frac{\mathfrak{F}}{\gamma }{A}_a[k_5{\overline{\mu }}_e(S,E_a)+k_7{\overline{r}}_{res}\left(E_a\right)]x_e$$

(18)

where $k_5$ and $k_7$ are constants.

3.8. Methane

Due to the low solubility of methane, it was assumed that this species passed directly to the gas phase of the anodic chamber once it was produced. The methane flow rate $q_{{\mathrm{CH}}_4}$ (mmol CH₄ d⁻¹) was:

$$q_{{\mathrm{CH}}_4}=k_2{A}_a{\overline{\mu }}_a\left(S\right)x_a$$

(19)

where $k_2$ is a constant.

3.9. Hydrogen

The species of hydrogen considered were the proton H⁺ and the molecular H₂. The protons were produced in the anodic chamber and traveled through the membrane to the cathode, where they were reduced:

$$2{\mathrm{H}}^{+}+2{\mathrm{e}}^{-}\to {\mathrm{H}}_2.$$

(20)

The hydrogen flow rate $q_{{\mathrm{H}}_2}$ (mmolH₂ d⁻¹) was:

$$q_{{\mathrm{H}}_2}={ϵ}_{cat}\frac{\gamma k_{{\mathrm{H}}_2}}{\mathfrak{F}}I$$

(21)

where ${ϵ}_{cat}$ (dimensionless) was the overall efficiency of the cathode for the hydrogen production and $k_{{\mathrm{H}}_2}$ = 0.5 (mmol H₂ mEQ⁻¹).

It was assumed that once the molecular hydrogen was produced, it passed immediately to the gas phase in the cathodic chamber due to its low solubility in water.

4. Simulations

The simulations of the model were carried out in MATLAB (2018a, MathWorks, Natick, MA, USA, 2018). The system presented in the Equations (14)–(17) is considered to have the next form:

$$\dot{\mathsf{\Phi }}=f(\mathsf{\Phi },u\left(t\right),\mathsf{\Pi })$$

(22)

where $\mathsf{\Phi }={[S,L_f,{\varphi }_a,{\varphi }_e]}^T$ is the vector of states; $f(\mathsf{\Phi },u(t),\mathsf{\Pi })$ is the vector field; $u\left(t\right)={[F\left(t\right),S^{in}\left(t\right),E_a\left(t\right)]}^T$ is the vector of inputs; and $\mathsf{\Pi }$ is the vector containing all the parameters of the model (such as those related to the kinetics, the thermodynamical constants, the yield coefficients, and the physical characteristics of the biofilm).

To integrate the system (22) in a determined time span, it is necessary to define the parameters of $\mathsf{\Pi }$ and the vector of initial conditions ${\mathsf{\Phi }|}_{t=0}$. The vector of inputs must also be defined, and it can be a function of time. Two types of simulations were done, chronoamperometry and voltammetry, as described and discussed in the following subsections.

4.1. Chronoamperometry

In Figure 1, Figure 2, Figure 3 and Figure 4 are shown the dynamical behavior of the system for 30 days of numerical simulation in which step change perturbations in $S^{in}$ and $E_a$ were carried out as shown therein. The parameters and operating conditions are shown in

Table 1. Set of nominal parameters.

Symbol	Value	Unit	Reference
$k_1$	$0.72$	mmol CH4 mmol S−1	[13]
$k_2$	$0.602$	mmol CO2 mmol S−1	[13]
$k_3$	$1.366$	mmol CO2 mmol S−1	estimate
$k_4$	$13.7$	mg VSe mmol S−1	[13]
$k_5$	$6.4$	mEQ mmol S−1	estimate
$k_6$	$1.393$	mmol CO2 mmol S−1	estimate
$k_7$	$0.177$	mEQ mg ${\mathrm{VS}}_{\mathrm{e}}^{-1}$	[17]
$k_8$	$0.044$	mmol CO2 mg VS−1	estimate
${\mu }_{a,max}$	$0.198$	mmol S mg ${\mathrm{VS}}_{\mathrm{a}}^{-1}$d−1	[15]
$K_{S,a}$	9.28 × 10−3	mmol S mL−1	[15]
$K_I$	0.256	mmol S mL−1	[15]
${\mu }_{e,max}$	0.132	mmol S mg VS−1d−1	[15]
$K_{S,e}$	3 × 10−5	mmol S mL−1	[15]
$E_{K_A}$	−0.156	V	[17]
$b_{in}$	$0.1$	d−1	[17]
$b_{res}$	$0.05$	d−1	[17]
$b_{det}$	60	cm−1 d−1	assumed
$D_e$	$0.753$	cm−2 d−1	[17]
$\gamma$	86,400	s d−1	estimated
$k_{bio}$	1 × 10−3	mA V−1 cm−1	[17]
${ϵ}_{cat}$	$0.9$	dimensionless	assumed
$k_{H_2}$	$0.5$	mmol H2 mEQ−1	assumed
$f_{exp}$	$0.8$	dimensionless	assumed
${\rho }_x$	1060	mg VS cm−3	[23]

and

Table 2. Set of additional parameters and operating conditions.

Symbol	Value	Unit
F	96,487	C mol ${\mathrm{e}}^{-1}$
R	8.314	atm L mol−1 k−1
P	1	atm
T	303	K
$V_a$	400	ML
F	400	mL d−1
$S_{in}$	$2000.0$	mg L−1
$A_a$	30.0	cm2
$L_f$	10	μm
$p_H$	7.0	dimensionless

. The initial operating condition for the simulation were: ${\mathsf{\Phi }}_0={[S,L_f,{\varphi }_a,{\varphi }_e]}^T$ = [2000 mg L⁻¹, 10 μm, 0.4, 0.3]^T. The following are some of the effects that can be observed in the simulations.

Figure 1a shows the dynamical behavior of the potential inside the biofilm and the potential applied ($z/L_f=0$) during the time simulation. The potential drop across the biofilm is significant because of the low conductivity value, even when the biofilm thickness was small. This drop was observed when the average potential (Figure 4c) was compared to $E_a$ ($z/L_f=0$). Figure 1b shows the dynamical behavior of substrate inside the biofilm during the time simulation. As can be seen, the concentration profile for each time step was practically the same along the biofilm. The numerical solution was computed based on the assumption that $\partial S^{bio}/\partial t=0$; thus, the distribution of $S^{bio}$ inside the biofilm immediately reached its steady-state. Furthermore, a high value of effective diffusion (based on previously-reported works) of the substrate in the biofilm ($D_e$) is proposed. Thus, as a result, similar values between local substrate concentration and substrate concentration in the bulk phase were obtained. Notice that the linear behavior of local potential along the biofilm remained only at the beginning of numerical simulation ($t<1$ d) concerning normalized biofilm thickness.

Figure 2 shows the local reaction rates during the time simulation. The local reaction rate ${\mu }_a\left(S^{bio}\right)$ (Figure 2a) behavior during $t<20$ d was considered constant. However, for $t>20$ d, the ${\mu }_a\left(S^{bio}\right)$ decreased exponentially with time. The local substrate concentration $S^{bio}$ did not vary along the biofilm (as mentioned before). Thus, as a result, the local reaction rate ${\mu }_a\left(S^{bio}\right)$ and the average reaction rate ${\mu }_a\left(S^{bio}\right)$ (Figure 3e) had a very similar dynamical behavior during the time simulation. The local reaction rates ${\mu }_e(S^{bio},E_a^{bio})$ and $r_{res}\left(E_a^{bio}\right)$ are a function of $E_a^{bio}$; and thus, as a result, these functions developed also a similar nonlinear dynamical behavior as $E_a^{bio}$ (Figure 1a). Notice that, the absolute error between the maximum or minimum local reaction rate and average reaction rate for ${\mu }_e(S^{bio},E_a^{bio})$ was bigger than $r_{res}\left(E^{bio}\right)$. However, the relative error was similar.

Figure 3 shows the average concentration of substrate inside the biofilm and the concentration in the bulk ($z/L_f=1$). Notice that both practically overlap. When the substrate concentration at the inlet of the system was perturbed, ($t=20$ d), the exoelectrogens microorganisms activity was not affected, i.e., the biomass concentration and the produced current remained constant. The primary factor here is that under these conditions, the multiplicative term of the growth kinetic related to the substrate was saturated. It can be proven that 95% of this term was reached when the acetate concentration was 19 mg L⁻¹. Therefore, the substrate concentration was not a controlling factor for the exoelectrogens microorganisms in the case of this study.

In this case, the kinetic growth rate of the exoelectrogen microorganisms was strongly affected by the potential. It was observed that in every moment, the average potential (Figure 4c) inside the biofilm was lower than $E_{K_A}=-156$ (mV). Hence, the Nernst–Monod term controlled the growth rate because this term was far from its saturation. It is important to mention that the effects on the biomass concentration (Figure 3c,d) and the current (Figure 4a) were evident when the applied potential changes (Figure 4c). If the potential increased, the exoelectrogens microorganisms transferred the electrons more easily, and their growth rate increased; thus, the biofilm thickness (Figure 3b) increased.

The produced electrons, according to the assumptions, came from two different sources, the metabolized acetate and the self-oxidized biomass, which was too small to be noticed. Moreover, not all the produced electrons went to the electrode. Indeed, some of them were used for the microorganisms’ growth and some others went to a different electron acceptor. The quantity $I_{ideal}$ (Figure 4a) was calculated under the assumption that all the produced electrons went to the electrode.

For the simulation, equal mass amounts of both microorganism types were assumed for the initial conditions. However, their concentration behaviors were completely different: the amount of exoelectrogens microorganisms increased, and that of the acetoclastic methanogens decreased (Figure 4d). Under the operational condition for $t>20$ d, there was no methanogen microorganism activity; thus, there was no methane production (Figure 4b). The main cause was the big difference between $k_1$ and $k_4$, parameters describing how much substrate is transformed into biomass. According to the assumed values, the exoelectrogens microorganisms produced almost 20 times more biomass per mole of acetate than the methanogens. It can also be observed that the biofilm thickness reached a steady-state, and this happened when the rate of detachment was equal to the biomass production. The inactivation process completely determined the steady-state amount of inert biomass in the biofilm.

4.2. Linear Voltammetry

Linear voltammetry is an electrochemical technique reported in many papers that tries to identify the behavior of the exoelectrogen microorganisms under the influence of an applied potential. In Figure 5 is presented a voltammetry simulation for different biofilm thicknesses. The process consisted of varying the applied potential and keeping the other inputs constant (Table 1 and Table 2). The scan was a rapid potential variation (in this case, 40 min) to assure that other changes in biomass concentration did not affect the overall behavior significantly during the voltammetry.

The simulated behavior was completely in line with previous experiments [26,30]. It is seen how the current increased for thicker biofilms (due to the more significant biomass concentration). In this simulation, the growth kinetics was affected only by the potential because of the term related to the substrate being saturated.

4.3. Additional Numerical Results

In the following, the case where the kinetics was not saturated is considered. It is important to notice that from a practical point of view, it is convenient to operate in saturation conditions to remove as much substrate as possible. However, in order to explore this nonsaturated scenario, two variables were modified. The first key variable was the inlet substrate concentration at the inlet of the system. In this sense, the inlet substrate concentration was reduced from $S_{in}=2000$ mg L⁻¹ to $S_{in}=800$ mg L⁻¹. The second modified parameter was the diffusion coefficient from $D_e=0.753$ cm⁻¹ d⁻¹ to $D_e=0.0753$ cm⁻¹ d⁻¹.

Figure 6a shows the dynamical behavior of potential inside the biofilm. Notice that the potential drop across the biofilm was significant, and notice as well that the average concentration in the biofilm was smaller than the bulk concentration. Figure 6a shows the dynamical behavior of substrate concentration inside the biofilm. In this case, the mass diffusion effect was significant at a low inlet concentration. For $t>15$ d, the concentration decreased from the bulk concentration of 200 mg L⁻¹ to zero.

Figure 7 shows the dynamical behavior of the reaction rates inside the biofilm. For $t>10$ d, as the mass diffusion effect was significant, the reactions rates varied significantly. In particular, for $t>10$ d ${\mu }_e$ (Figure 7b) and $r_{res}$ (Figure 7c), they were near to zero from $z/L_f=1$ to $z/L_f=0.5$. Figure 8a, shows that there was a significant difference between bulk substrate concentration and substrate average along the biofilm. Under this scenario, the numerical results of biomass fraction, thickness grown, and biomass concentration (see Figure 8b–d, respectively) were lower than those obtained using a typical set of parameters obtained from the literature. Such a difference was significant for ${\overline{\mu }}_a$ (see Figure 8e). As a result, the lower current was obtained (see Figure 9a).

At a low diffusion coefficient and low substrate concentration at the inlet of the system, the profile of the substrate along the biofilm (Figure 6b) changed concerning the parameters obtained from the literature (Figure 1b).

5. Perspectives and Recommendations

MEC process has been proven to be a promising technology with the potential for a plethora of applications [8,31,32]. Indeed, existing MEC technologies can be further refined, and new applications can arise [31,33]. For instance, the theoretical framework concerning the proposed biofilm model may be taken as a basis for other applications like chemicals’ synthesis or treatment of recalcitrant compounds. However, the scaling up of MEC systems from lab scale to pilot scale and larger scales is still an open challenge that requires a multidisciplinary approach [32,34]. Improving the modeling and understanding of MEC systems is crucial to establish an optimal operational condition, i.e., closed loop operation [32]. The proposed dynamical model constitutes not only a suitable, but also a promising alternative to the actual necessity of the mathematical description of MEC systems for optimization and control purposes, in a process alone or a technological building block within the concept of the environmental biorefinery [35].

6. Conclusions

In this study, a model was proposed and simulated with parameters from the literature. The results of the simulations were in line with the reported behavior of this system in previous papers.

Some of the most important results that can be drawn from this model and its simulations are:

• The biofilm reached a steady-state thickness [22,27]
• Even when the biofilm thickness was small, the charge transfer through the film played a significant role because of the potential losses [26,30]
• The response of the system to different perturbations in the applied potential and the inlet substrate concentration accorded with the experimental data of [36]
• It was observed that the potential loss and the biofilm processes were the controlling factors.