A 2D Modelling Approach for Predicting the Response of a Two-Chamber Microbial Fuel Cell to Substrate Concentration and Electrolyte Conductivity Changes

Abstract

Bioelectrochemical systems have been the focus of extensive research due to their unique advantages of converting the chemical energy stored in waste to electricity. To acquire a better understanding and optimize these systems, modelling has been employed. A 2D microbial fuel cell (MFC) model was developed using the finite element software Comsol Multiphysics^® (version 5.2), simulating a two-chamber MFC operating in batch mode. By solving mass and charge balance equations along with Monod–Butler–Volmer kinetics, the operation of the MFC was simulated. The model accurately describes voltage output and substrate consumption in the MFC. The computational results were compared with experimental data, thus validating the model. The voltage output and substrate consumption originating from the model were in agreement with the experimental data for two different cases (100 Ω, 1000 Ω external resistances). A polarization curve was extracted from the model by shifting the external resistance gradually, calculating a similar maximum power (47 mW/m²) to the observed experimental one (49 mW/m²). The validated model was used to predict the MFC response to varying initial substrate concentrations (0.125–4 g COD/L) and electrolyte conductivity (0.04–100 S/m) in order to determine the optimum operating conditions.

1. Introduction

Bioelectrochemical systems (BES) have been the focus of research due to their ability to convert the chemical energy contained in various wastewaters to electricity [1,2]. The operation of these systems is based on the oxidation of an electron donor followed by the reduction of an electron acceptor. The oxidation reaction takes place in the anode chamber and is catalyzed by electrogenic bacteria, which have the ability to transfer the electrons to their exterior [3]. Similarly, the reduction reaction occurs in the cathode chamber and may be catalyzed by microorganisms (biocathode) or catalysts such as platinum, activated carbon, etc. (abiotic cathodes) [4]. Different types of BES have been developed, such as microbial fuel cells (MFCs), microbial electrolysis cells (MECs), microbial solar cells (MSCs) [5], microbial electrosynthesis cells [6], microbial desalination cells (MDCs) [7], and enzymatic fuel cells (EFCs) [8]. These electrochemical cells have been developed for simultaneous wastewater treatment, current generation, and in the case of MECs production of other substances such as hydrogen [9]. The most studied BES systems are Microbial Electrolysis Cells (MECs) and Microbial Fuel Cells (MFCs). In the case of MEC systems, an external voltage is applied in order for electrolysis to take place. [10]. On the other hand, the redox reactions that take place in a microbial fuel cell (MFC) occur spontaneously, generating electricity by oxidizing an organic substrate [11]. The plethora of bioelectrochemical processes which are catalyzed by microorganisms are exhibited in the multiple BESs that have been developed.

The phenomena that take place in BESs are varied and complex. The most important phenomena that take place in BESs have been presented and summarized by Rabaey et al. [12]. The different processes that occur within BESs, such as MFCs, can be broken down into the following phenomena: (i) microbial metabolism on the anode electrode; (ii) diffusion of the organic substrate from the bulk into the biofilm of the anode chamber initiating the oxidation reaction of the electron donor; (iii) transport of the organic substrate by diffusion and by stirring (convection) in the bulk of the anolyte; (iv) oxidation reaction and production of protons and electrons; and (v) transfer of electrons through the wiring connection (external resistance) from the anodic electrode to the cathodic electrode. In the case of an MEC, instead of the external resistance a voltage is applied between the anode and the cathode. The protons are diffused into the anode and then transported via convection and diffusion through the separator to the cathode chamber (vi), to either participate in the reduction reaction or to balance the concentration gradient between the anode and the cathode. Other substances produced by oxidation are diffused into the bulk of the anolyte. In the cathode chamber, the electron acceptor comes into contact with the cathodic electrode and, combined with the electrons, starts the reduction reaction (vii). The products of the reduction are deposited on the cathode electrode, settled in the cathode chamber, or mixed with the bulk of the catholyte (viii). Moreover, in the case of a bio-cathode, microbial metabolism can be studied in the cathode. Overall, charge transfer and balance are examined in the BES, along with losses regarding the developed voltage [12].

Modelling is a powerful tool in seeking to understand the complex processes occurring within BES systems. Several models have been proposed in this direction. Because of the various MFC setups (single chamber, dual chamber) and the variety of materials used for the electrodes and the catalysts, in combination with the different phenomena during the MFC operation, many models have been developed [13,14]. One of the first models describing MFC operation was created by Zhang and Halme in 1995 [15]. Zhang and Halme [15] simulated the operation of a dual-chamber MFC using Monod kinetics to describe the bacterial growth and substrate consumption, coupling this with Faraday’s law to calculate the current generation. Their model was based on the assumption that mass transport and the cathodic reaction are faster than the biochemical and oxidation reactions, and thus focused on the latter. As one of the first attempts at MFC modelling, they highlighted the importance of maintaining control over the process through simulation [15]. An example of a model combining bioelectrochemical kinetics along with mass and charge balances in a dual-chamber MFC was developed by Zeng et al. in 2009 [16]. The developed model aimed to describe an acetate-fed dual chamber MFC with the ability to adapt to other organic substrates. In order to describe the electrochemical phenomena, a combination of the Butler–Volmer expression and Monod kinetics was utilized. Several of the parameter values were extracted from experimental data, while others were estimated by the model. Zeng et al. concluded that the cathodic reaction is the liming factor for current generation in an MFC. Furthermore, their results indicated that an increase in the electron donor concentration was effective in boosting MFC performance [16].

Various approaches have been presented in the literature, with different focus points. A study on the biofilm of electrochemically active bacteria during the operation of an MFC was conducted by Belleville et al., 2019 [17]. A 2D model was developed, examining at a cellular level the microorganisms on the anodic electrode. Taking into account microbial growth and segregation, their model aimed to predict the different types of bacteria (methanogens, electrogens) that grow, perish or detach from the electrode during the operation of a glucose-fed MFC [17]. Serra et al., 2019 focused on the power output of MFCs through polarization curves and created a steady-state electrical model to simulate them [18]. The polarization experiments were conducted by varying the external resistance on six different MFCs from 1000 Ω to 20 Ω, with simultaneous refreshment of the synthetic wastewater. The model had different kinetics for the different power losses (activation, ohmic, concentration) which may be present during MFC operation [18].

Furthermore, the various operation modes of MFCs have been examined. A 2D model for simulating the operation of a continuously fed, single-chamber MFC was developed by Day et al., 2021 [19]. Different parameters such as the variation of the hydraulic retention time (HRT) and the alteration of MFC geometry were examined in order to determine the optimal conditions for simultaneous maximum current generation and substrate consumption. Matsena et al., 2021 created a model simulating the operation of an MFC with hexavalent chromium as the electron acceptor [20]. Monod kinetics were used to describe the growth of the biofilm taking into account inhibition by the substrate and the intracellular mediator. Moreover, two types of microorganisms were considered: one that contributes to electricity generation and one that consumes the substrate without releasing electrons. To describe the electrochemistry, the Butler–Volmer expression was utilized, while Faraday’s law was used to couple the electrochemistry with the microbial activity. In order to simulate the reduction of hexavalent chromium, they combined Faraday’s law with Butler–Volmer, and the mass balance equation was derived from the reaction of dichromate (Cr₂${\mathrm{O}}_7^{-}$) to Cr(III) [20].

Apart from simulations describing MFC operation, models have been developed to support experiments. Thus, Sindhuja et al., 2016 conducted electrochemical impedance spectroscopy (EIS) experiments in a dual chamber MFC and developed a model to extract the equivalent circuit [21]. The model was developed for an MFC with charcoal electrodes, and was used along with a fitting model of a Nyquist plot for graphite electrodes. Their aim was the determination of internal resistances in both experimental configurations (charcoal, graphite based electrodes respectively) and the estimation of the rate determining step in their MFC operation with glucose as the organic substrate [21]. Oliveira et al., 2013 published a study presenting a 1D MFC model [22]. This particular model incorporated heat, charge, and mass transfer phenomena at a steady state across an MFC configuration consisting of an anodic electrode, a biofilm, the anolyte solution, the separator (proton exchange membrane), the catholyte, and the cathodic electrode. They successfully compared their computational data to the respective results presented by Zeng et al., 2009 [16,22]. By briefly presenting these models, it is easy to deduce that MFC modelling is applied over different focus points (i.e., cellular level, spatial and time dependent, or dimensionless), resulting in a plethora of publications describing different aspects of this technology.

The aim of this work was to develop a time-dependent 2D MFC model containing mass and charge conservation and transfer phenomena along with combined Butler–Volmer–Monod kinetics. The model was developed with minimal computational requirements in order to extract quick and accurate results taking into account the geometry of the cell and the materials that comprise it. This model examines glucose consumption in an H-type MFC, calculating the voltage output, substrate consumption and polarization curves. The model proposed an approach to quickly predict and validate optimal conditions for MFC operation. This was accomplished by the use of electrochemical kinetics along with parameters determined by experimental data fitting. The focus of this work was the creation of an improved model to quickly predict accurate conditions, in order to improve the performance of the particular MFC setup.

2. Materials and Methods

2.1. Experimental Setup

The model was based upon a dual-chamber MFC in order to validate the computational results with experimental data. Specifically, an H-type MFC was operated in batch mode, consisting of an anode (0.3 L) and a cathode (0.3 L). The solutions in both chambers were continuously stirred. Plain graphite paper (3.8 cm × 2.5 cm) and graphite cloth (3.8 cm × 2.5 cm) coated with Pt as the oxygen reduction catalyst were used as the anodic electrode and the cathodic electrode, respectively. Titanium wire and a silver conductive epoxy (Conductive Epoxy, Circuit Works) were used for the electrode connection. A proton exchange membrane (Nafion^® 117, DuPont–PEM) was placed between the two chambers. Both chambers were under continuous stirring and operated in a temperature control environment set at 30 °C. A resistance box set at 100 Ω was connected to the cell. The cell was operated with synthetic glucose wastewater (1 g COD/L) in the anode, containing a phosphate buffer (3.67 g/L NaH₂PO₄ and 3.45 g/L Na₂HPO₄), potassium chloride (0.16 g/L KCl), sodium bicarbonate (5 g/L NaHCO₃), and solutions containing elements necessary for the microbia (1% v/v, presented in

Table 1. Composition of trace elements solutions used in the synthetic glucose feed.

Component	Concentration (mg/L)
Solution A
CaCl2·2H2O	22,500
NH4Cl	35,900
MgCl2·2H2O	16,200
KCl	117,000
MnCl2·4H2O	1800
CoCl2·6H2O	2700
H3BO3	513
CuCl2·2H2O	243
Na2MoO4·2H2O	230
ZnCl2	189
NiCl2·6H2O	200
H2WO4	10
Solution B
FeSO4	700
Solution C
(NH4)2PO4	7210

, based on [23]). During the acclimation period, anaerobic sludge (10% v/v) was added to the synthetic glucose wastewater (10 g COD/l). For the cathode phosphate buffer, (3.67 g/L NaH₂PO₄ and 3.45 g/L Na₂HPO₄) and potassium chloride (0.16 g/L KCl) were used. The cathode chamber was continuously sparged with air, supplying the electron acceptor (O₂).

The cell’s potential was recorded at set intervals (2 min) using an Agilent Keysight 34972A LXI Data Acquisition/Switch Unit. The pH and conductivity were measured by digital instruments, (WTW INOLAB PH720) and (WTW INOLAB), respectively. Soluble COD was measured according to standard methods [24]. Polarization experiments were conducted by altering the resistance on the resistance box. Initially, the system achieved open circuit voltage (OCV) by removing the external resistance for 3 h, thus resulting in infinite resistance between the two electrodes. Then, the resistance box was connected to the cell and different consecutive resistances were applied (1 MΩ–0 Ω). Each external resistance was applied for 6 min, recording the voltage and the current every 2 min. The current was measured using a multimeter.

2.2. Model Description

The model was developed in the finite element software (FEM) Comsol Multiphysics^® Version 5.2. The model focused on the description of the voltage output, the organic substrate (glucose) consumption, and the estimation of the maximum power through polarization curves. Experimentally, the glucose concentration was measured in g COD/L; however, in the model, the corresponding glucose concentration (mol/m³) was determined. For this reason, the organic substrate results were normalized with the initial value (C/C₀).

The geometry was a vertical section (2D) of the MFC, depicted in Figure 1a. It consisted of an anode, a cathode, the separator (PEM membrane [25]), and the two electrodes (2.5 cm × 3.8 cm). The third dimension (z-axis) was considered equal to the width of the electrodes (2.5 cm). The distance between the two electrodes was 15 cm. The cell was considered to operate in batch mode and at a set temperature (32 °C). As the electron donor to be oxidized by the electrogen biofilm, glucose was chosen. The reaction assumed to take place in the biofilm was [26]:

C₆H₁₂O₆ + 6 H₂O → 6 CO₂ + 24 H⁺ + 24 e⁻

(1)

The biofilm was considered to be on the surface of the anodic electrode such that glucose reacted as it came into contact with it. Despite the porosity of graphite paper and the presence of microorganisms, the reaction was assumed to take place homogeneously on the electrode. In the anode chamber, the carbon dioxide was considered dissolved in the liquid, thus taking into account only the liquid phase. Furthermore, due to the presence of a phosphate buffer, the changes in the anodic pH were considered negligible. The concentration of microorganisms was considered constant, as the assumption was made that the growth rate was equal to the decline rate. This assumption led to a constant biofilm present on the anodic electrode. Moreover, the biofilm was assumed to consist only of electrogenic bacteria and no competing biomass was considered to be present. The description of the electrochemical reaction rate (mol/m²/s) taking place in the biofilm used the Monod–Butler–Volmer expression. This equation combines the organic substrate consumption rate and the effect of electrochemical phenomena on it, taking into account the overpotential:

$$R_1=k_1\frac{C_{Glucose}}{K_{Glucose}+C_{Glucose}}\exp \left(\frac{a_{a}F}{RT}{\eta }_a\right)$$

(2)

where k₁ (mol/(m²·h)) is the maximum specific growth rate multiplied by the biomass concentration (k₁ = $k_1^0·X_{bio}$), C_Glucose (mol/m³) is the concentration of glucose in the anode chamber, K_Glucose (mol/m³) is the half velocity rate constant for glucose, $a_a$ is the anodic transfer coefficient, F (C/mol) is Faraday’s constant, R (J/(mol·K)) is the gas constant, T (K) is the temperature, and ${\eta }_a$ (V) is the anode overpotential. No inhibition of microbial activity by the organic substrate was considered.

Through the separator, only protons were transferred. The cathodic chamber was considered to be operating with a phosphate buffer as the catholyte, resulting in negligible variations of pH. The electrons were transferred to the cathode electrode through external resistance, which was set at 100 Ω during batch cycles and varied accordingly during polarization experiments. The dissolved oxygen was assumed to react when in contact with the cathode electrode. The concentration of dissolved oxygen was considered constant because of the continuous supply of air in the cathode chamber. The reduction reaction (proton consuming [12]) taking place on the cathodic electrode was

O₂ + 4 H⁺ + 4 e⁻ → 2 H₂O

(3)

It was assumed that no diffusion of oxygen from the cathode to the anode would take place. To describe the rate of the reaction in Equation (3), a Monod–Butler–Volmer kinetic was employed, assuming a Monod-type dependence on dissolved oxygen concentration [27]:

$$R_2=k_2\frac{C_{Oxygen}}{K_{Oxygen}+C_{Oxygen}}\exp \left(\frac{a_{C}F}{RT}{\eta }_c\right)$$

(4)

where k₂ (m¹²/(m⁴·h)) is the forward rate constant of cathode reaction, C_Oxygen (mol/m³) is the concentration of dissolved oxygen in the cathode chamber, K_Oxygen (mol/m³) is the half velocity rate constant for oxygen, $a_C$ is the cathode transfer coefficient, and ${\eta }_c$ (V) is the cathode overpotential.

The correlation between reaction and rate and current development was carried out using Faraday’s law:

$$R_i=\frac{{\nu }_i{i}_{A/C}}{n_{i}F}$$

(5)

where ν_i is the stoichiometric coefficient (i = glucose or oxygen), i_A (A/m²) is the current density (anode or cathode), and n_i the number of electrons that take part in the reaction.

The algorithm of the model is presented in Figure 1b. After selection of an initial glucose concentration, the simulation is initiated and the reaction rate on the anode electrode surface is calculated by Monod–Butler–Volmer kinetics. The initial value for the voltage between the two electrodes is 0 V. Using Faraday’s law, the current density on the anode is calculated based on the reaction rate. Through external resistance, the current brings electrons to the cathode electrode, where oxygen reduction takes place. The respective current density is calculated by Faraday’s law based on the oxygen reduction reaction rate. The charge balance equations calculate the total transfer of charge throughout the cell and the voltage developed between the two electrodes. The overpotential is determined based on the voltage and the standard potentials of the reactions. Simultaneously, the decrease in the glucose concentration and the respective glucose distribution is determined. As a new glucose concentration on the electrode surface is calculated, a new reaction rate is determined, taking into account the new overpotential value.

2.3. Mass and Charge Transfer Equations

The anode and the cathode chambers were considered to be full of the aforementioned synthetic solutions. For the transfer of mass through the anode and the cathode the Nernst–Planck equation was used. This equation takes into account the change of concentration in time and the transport of chemical species by diffusion, convection, and migration. This equation takes into account the transport of species due to the concentration gradient (diffusion), by the movement of the bulk of the fluid (convection), and due to the presence of an electric field (migration).

$$\frac{\partial \mathrm{c}}{\partial \mathrm{t}}-\nabla ·\left[D\nabla \mathrm{c}-u\mathrm{c}+\frac{zF}{RT}c\left(\nabla \phi \right)\right]=\{\begin{array}{c}0, in the bulk \\ R, on the electrode surface\end{array}$$

(6)

where D (m²/s) is the diffusion coefficient of the species, u (m/s) is the velocity vector of the liquid, z is the valence of the ionic species, φ (V) is the potential, and R (mol/(m²·h)) is the reaction rate. A no flux condition was applied on the walls; similarly, no flux was applied on the separator. Inside the cell, an incompressible Newtonian fluid was assumed.

For the charge transfer in the cell, Ohm’s law was used along with charge balance equations. More specifically, a uniform electrolyte medium was assumed for the transfer of charged ions, and activation overpotentials were taken into account as well:

$$\nabla ·i_i=Q_i$$

(7)

$$i_i={\sigma }_i\nabla {\phi }_i$$

(8)

where i_i (A/m²) is the current density (i = electrode/electrolyte), Q (C) is the total charge, σ (S/m) is the conductivity, and φ (V) is the potential. Ohmic changes in the electrical connections were considered negligible. For the calculation of the cell’s voltage, the standard reaction potential, the anode and cathode overpotential, and the ohmic losses in the electrolyte were considered using the following equation:

$${\mathrm{V}}_{cell}=E^0+{\eta }_{anode}+{\eta }_{cathode}+R_{cell}{i}_{cell}$$

(9)

where E⁰ (V) is the standard redox potentials and R_cell (Ω) is the internal resistance of the cell, consisting of the separator and the electrolyte resistances. Furthermore, Ohm’s law was used to calculate the voltage between the anode and cathode electrodes:

$${\mathrm{V}}_{cell}=I{R}_{ext}$$

(10)

where I (A) is the current of the closed circuit and R_ext (Ω) is the external resistance connected in the MFC. After the initial simulation runs, the capacitance of the electrodes was incorporated into the model in order to further examine its effect. The following equation was added to take into account the changes in the potential in the electrode electrolyte interface, as well as the capacitance:

$$i=\left(\frac{\partial \left({\phi }_{electrode}\right)}{\partial t}\right)C_{electrode}$$

(11)

where C_electrode (F/m) is the capacitance of the electrode. In the anode, the capacitance takes into account the biofilm’s existence and is different for every acclimation, namely every biofilm. The capacitance depends on the type of the electrode as well.

The parameters of the aforementioned equations were extracted by a fitting model developed using Aquasim 2.0 software [28]. This model was a modified version of the one used by Zeng et al., 2009 [16]. The model contained altered mass balance equations in order to simulate a batch reactor instead of a CSTR. Glucose was assumed as the electron donor, and Andrews kinetics were used instead of Monod. The resulting Andrews–Butler–Volmer equation is as follows:

$$R_i=k_i\frac{C_{Glucose}}{K_{Glucose}+C_{Glucose}+\frac{C_{Glucose}{}^2}{K_i}}\exp \left(\frac{a_{a}F}{RT}{\eta }_a\right)$$

(12)

where K_i (mol/m³) is the inhibition of the glucose concentration. Experimental data on the operation of the previously-described dual chamber MFC were used as input values for the estimation of the parameters. The results and additional details of the Aquasim model are presented elsewhere [28]. The values of the parameters used are presented in Table 2.

3. Results

3.1. Model Results and Validation

The model presented fast convergence in under one minute, simulating the consumption of glucose in the MFC and the simultaneous voltage output. Initially, the results from the model were compared with the operation of the two-chamber MFC, with the external resistance set at 100 Ω. Figure 2a presents the concentration of the organic substrate calculated during one batch cycle as extracted by the model, as well as the concentration obtained from the experiments. Figure 2b shows the voltage output of the cells as calculated by the model along with the respective values originating from the experiment. Specifically, the concentration was normalized to the initial substrate concentration (1 g COD/L), while voltage recording started as soon as the MFC was fed with fresh glucose synthetic wastewater.

As shown in Figure 2a, the concentration initially presents a linear drop, then tends asymptotically to zero (blue line Figure 2a). The COD removal achieved by the cell was 91% (Figure 2a). The simulation continued until the substrate concentration inside the cell reached 0 mol/m³. The operation cycle lasted for 75 h, while the simulated time for the model was 100 h (Figure 2a). The model predicted a 98% substrate consumption at the 75 h mark. The experimental measurements of substrate concentration (red circles) were in agreement with the simulation data.

The voltage output originating from the model at 100 Ω was close to the experimental data (Figure 2b). In particular, the maximum voltage was 29 mV, while the respective voltage peak during the experiment was 31 mV. The voltage plateau was close, and with the same duration for both experimental and computational results (~40 h). The difference between the experiment and the model was that the plateau formed immediately in the model (0 h), while in the experiment the maximum voltage was achieved after 10 h of cell operation. This was attributed to glucose diffusion in the biofilm and the activation energy required to initiate glucose oxidation. On the other hand, in the model the reaction takes place on the electrode surface and comes immediately into contact with the organic substrate, thus achieving the maximum value at the beginning of the cycle. The second difference between the results of the experimental data and the model prediction is observed after the plateau and the voltage decrease trend. The experimental data show that it takes 19 h for the voltage to be reduced to 0 V; however, the model requires 59 h to decrease this value to 0 V. This deviation is attributed to the ohmic losses which are present during the experiment because of the electrical connections. In the case of the model (100 Ω), a slower and smoother voltage decline takes place. Overall, the model fitting regarding the voltage output of the units at 100 Ω external load is considered to be satisfactory. The experimental data used in the modified Zang model {Formatting Citation} were extracted from the operation of the cell with 100 Ω external resistance. The voltage output and the substrate consumption were in good agreement with the experimental measurements in the MFC operated with 100 Ω, as these experimental data were used for the parameter calculation.

For better understanding of the organic substrate distribution, four different images are presented in Figure 3. The concentration of glucose was normalized to the initial value (C₀).

The distribution of the substrate is presented in Figure 3a, at a point when the oxidation had not started yet. At 25 h (Figure 3b) the substrate consumption was 45%, at 50 h (Figure 3c) it was 87%, and at 75 h it was 98% (Figure 3d). The distribution of the substrate in the anodic chamber is uniform, with the exception of a gradient present at 25 h and 50 h. This was observed around the anodic electrode, as expected, as the reaction takes place on its surface. Moreover, in all four cases (Figure 3a–d) the effect of the separator was seen to inhibit the transfer of glucose from the anode to the cathode (C_cathode = 0, dark blue color).

In order to validate the model and test its ability to predict MFC performance with R_ext 1000 Ω, more data were compared. Figure 4 presents the results of the simulation of glucose concentration and the voltage output during cell operation (at R_ext 1000 Ω) in comparison with the respective experimental data.

By shifting the external resistance in the model to 1000 Ω, similar results were obtained with the respective experimental data (Figure 4). The experimental COD removal was 89% at the end of the batch cycle. The cycle duration for the 1000 Ω experiment was 92 h, close to the simulated time (110 h) for the full depletion of substrate in the anode chamber. At 92 h, the model calculated 98% substrate consumption (Figure 4a). A similar pattern is observed from the results originating from the system set at 100 Ω, as expected, as the consumption of glucose primarily depends on the Monod kinetics and not on the external resistance. The substrate concentration decreased similarly in the model and the experiment. The voltage output peak was 220 mV for the model and 210 mV for the experiment (Figure 4b). The voltage plateau was maintained for 63 h in the simulation and 72 h in the experiment. The voltage output reached its maximum value (210 mV) after 10 h of cell operation, while the model predicted the voltage to reach its maximum value (220 mV) at 0 h. The voltage decrease lasted for 52 h simulated time, while this decrease in the experiment was faster (7 h). Similar to Figure 2b, the time difference between the model and the experiment is attributed to glucose diffusion within the biofilm, and to the voltage losses on account of the electrical connections and cell geometry (H-type). Moreover, when the simulated voltage plateau ceased in both cases, at 40 h and 63 h, respectively (see black arrows in Figure 2b and Figure 4b), the organic substrate was reduced by 75% and 81%, respectively (see black arrows in Figure 2a and Figure 4a). The prediction of the MFC operation with an external resistance set at 1000 Ω was considered successful, validating the developed model.

After simulating the MFC operation with two different external resistances (100 Ω, 1000 Ω) and comparing the results with experimental data, the next step was the extraction of the polarization curve. A parametric study was carried out on the external resistance variable, testing a range of values (0.5 MΩ–0 Ω).

The polarization curve calculated by the model (blue line Figure 5) was in agreement with the respective curve resulting from the experiment. Both power and current values were normalized to the surface of the anodic electrode (Table 2) in order for the results to be easily comparable with other cases. The model’s maximum power density (47 mW/m²) was obtained at an external resistance equal to 5000 Ω and current density equal to 63 mA/m². The corresponding experimental power density (49 mW/m²) was achieved at 1900 Ω and 80 mA/m². Despite the apparent consistency of these results, the same maximum power was achieved for different external resistances (Simulation: 5000 Ω, Experiment: 1900 Ω) and consequently different current densities (see blue and red arrows in Figure 5). On the other hand, the model was able to predict the range of the current densities 2.6–154 mA/m² produced by the MFC (0–149 mA/m²). The difference between the experimental and the computational polarization curves was observed in the 20–140 mA/m² range. In Figure 5, the I–V curves for both the experiment (red line) and the simulation (blue line) are presented. The slope of the lines indicates in both cases the ohmic resistance as the main cause of electrochemical losses. The OCV of the model was 0.75 V, while the respective value measured in the experiment was 0.7. In the simulated I–V curve, apart from the ohmic losses the effect of substrate diffusion is visible for high current densities and low voltages. The experimental I–V line does not indicate other electrochemical losses apart from ohmic.

Table 2. List of model parameters presented in the order their equations were presented (1–12).

Symbol	Description	Unit	Value	Source
$k_1$	Maximum specific growth rate	mol m−2 h−1	6 × 10−3	Calculated by Aquasim model [28]
$K_{Glucose}$	Half–velocity rate constant for glucose	mol m−3	3 × 10−4
$a_a$	Anode transfer coefficient	–	0.05
$k_2$	Forward rate constant of cathode reaction	m12 mol−4 h−1	9.19 × 10−5
$K_{Oxygen}$	Half–velocity rate constant for oxygen	mol m−3	4 × 10−3
$a_C$	Cathode transfer coefficient	–	0.7
$C_{anode}$	Anode Capacitance	F m−2	13,721
$C_{cathode}$	Cathode Capacitance	F m−2	500
$K_S$	Glucose inhibition constant	mol m−3	37 × 10−3
F	Faraday’s constant	Coulombs mol−1	96,485	[29]
R	Gas constant	J mol−1 K−1	8.31
$C_{Oxygen}$	Dissolved oxygen concentration in the cathode chamber	mol m−3	0.3125	[30]
$D$	Glucose diffusion coefficient	m2 s−1	0.5 × 10−9
${\sigma }_i$	Electrolyte conductivity	S m−1	1.2	Experimental values
${\sigma }_i$	Electrode conductivity	S m−1	10
$E^0$	Open Circuit Voltage	V	0.75
$S_{electrode}$	Electrode surface	cm2	19
$S_{separator}$	Separator surface	cm2	1
$d_{electrode}$	Electrode distance	cm	17

Table 2. List of model parameters presented in the order their equations were presented (1–12).

Symbol	Description	Unit	Value	Source
$k_1$	Maximum specific growth rate	mol m−2 h−1	6 × 10−3	Calculated by Aquasim model [28]
$K_{Glucose}$	Half–velocity rate constant for glucose	mol m−3	3 × 10−4
$a_a$	Anode transfer coefficient	–	0.05
$k_2$	Forward rate constant of cathode reaction	m12 mol−4 h−1	9.19 × 10−5
$K_{Oxygen}$	Half–velocity rate constant for oxygen	mol m−3	4 × 10−3
$a_C$	Cathode transfer coefficient	–	0.7
$C_{anode}$	Anode Capacitance	F m−2	13,721
$C_{cathode}$	Cathode Capacitance	F m−2	500
$K_S$	Glucose inhibition constant	mol m−3	37 × 10−3
F	Faraday’s constant	Coulombs mol−1	96,485	[29]
R	Gas constant	J mol−1 K−1	8.31
$C_{Oxygen}$	Dissolved oxygen concentration in the cathode chamber	mol m−3	0.3125	[30]
$D$	Glucose diffusion coefficient	m2 s−1	0.5 × 10−9
${\sigma }_i$	Electrolyte conductivity	S m−1	1.2	Experimental values
${\sigma }_i$	Electrode conductivity	S m−1	10
$E^0$	Open Circuit Voltage	V	0.75
$S_{electrode}$	Electrode surface	cm2	19
$S_{separator}$	Separator surface	cm2	1
$d_{electrode}$	Electrode distance	cm	17

3.2. Different Initial Substrate Concentrations

The model was further validated in terms of its ability to predict MFC performance when the initial glucose concentration is changed. The initial glucose concentration ranged between 0.125–4 g COD/L. In terms of Chemical Oxygen Demand removal (COD), 90% COD removal was considered adequate in a well-performing MFC. The time required for the 90% COD removal was defined as the time required for the MFC model to achieve a satisfying substrate treatment. Furthermore, the duration of the maximum voltage output was maintained for the different initial substrate concentrations. The results are presented in Figure 6.

The duration of the plateau during which maximum voltage output was maintained increased with increasing glucose concentration (Figure 6, black line). A similar increase was observed for the time to 90% COD removal, corresponding to a glucose increase (Figure 6, red line). In the case with the glucose concentration at 0.125 g COD/L the maximum voltage value was 15 mV, which was lower than that achieved for the other concentrations (29 mV), and was maintained for less than 1 h. For 0.25, 0.3 g COD/L the maximum voltage (29 mV) was maintained for 1 h as well. As the initial glucose concentration increased, the duration of the maximum voltage increased as well, presenting a linear correlation between the two. For the 90% COD removal duration, a similar pattern was extracted from the model. The time required for 90% COD consumption was higher than the respective duration of the maximum voltage required for lower initial concentrations (0.125–2 g COD/L). On the other hand, for values higher than 2 g COD/L the two lines intersected at approximately 2.5 g COD/L and the duration of the maximum voltage was higher than the corresponding time required for 90% COD removal. This result indicates that there may be an initial concentration which, when fed to an MFC will achieve a more intensive and effective operation in successive batch cycles. Moreover, for lower initial concentrations the MFC underperforms in terms of voltage output, as it did not achieve the maximum voltage value. The maximum voltage (29 mV) achieved did not depend on the increasing initial glucose concentration for the larger part of the range examined (0.4–4 g COD/L).

3.3. Different Initial Electrolyte Conductivities

In order to examine the model’s capabilities, a parametric study was conducted on the effect of electrolyte conductivity, maintaining all other parameters at their respective initial values (Table 2) and the external resistance at 100 Ω. The electrolyte conductivity was adjusted in the experiments by the addition of potassium chloride, sodium hydroxide, and trace elements, which altogether increased the conductivity of the synthetic glucose solution to 1.2 S/m. The range of the electrolyte values tested was 0.036 S/m–100 S/m. The initial glucose concentration corresponded to 1 g COD/L.

Figure 7 presents the results of the electrolyte conductivity parametric study. The maximum voltage output was 30 mV, and the lowest voltage was 2.1 mV. For a better presentation of the computational results, a logarithmic scale was used for the x-axis, aiming to indicate the effect of the electrolyte conductivity on the voltage output. The point of the electrolyte value (1.2 S/m) used in the experiments and the initial runs (100 Ω, 1000 Ω) is highlighted in Figure 7 (black arrow). An increase in the maximum voltage output was observed as electrolyte conductivity increased. Furthermore, at higher conductivity values (>2 S/m) the maximum voltage from each simulation converged to the same value (30 mV). The initially-selected value of the electrolyte conductivity (1.2 S/m) achieved a similar maximum voltage (29 mV). These results indicate that by increasing the electrolyte conductivity up to 2 S/m, the maximum voltage output of the MFC increased as well. Further increasing the electrolyte conductivity (>2 S/m) did not have a similar effect on voltage, and the maximum value remained the same. Based on this conclusion, it was deduced that in order to increase the performance of an MFC the addition of electrolyte salts was effective up to a critical value, beyond which adding extra conductivity boosters had no effect on the maximum voltage.

The different electrolyte conductivities did not have an effect on COD removal. This observation was expected as well, as conductivity does not directly affect the substrate consumption, rather affecting the voltage which, through the overpotential, is implemented in Monod–Butler–Volmer. The voltage was not high enough to have an impact on glucose consumption, nor did it affect the transport of species due to migration, as the electric field was weak.

4. Conclusions

In this work, a 2D MFC model was developed in order to simulate the complex operation of MFCs. Electrochemical kinetics along with mass and charge transfer equations were solved in the FEM software Comsol Multiphysics^®. The results were compared with experimental data from the operation of the MFC on which the model was based. The model effectively predicted the value of the maximum voltage output of MFCs operated at R_ext 100 Ω (simulation: 29 mV; experiment: 30 mV) and R_ext 1000 Ω (simulation: 220 mV; experiment: 210 mV). Moreover, by shifting the external resistance of the MFC, the polarization curve extracted from the model resulted in a similar maximum power (47 mW/m²) as the corresponding experimental data (49 mW/m²). By changing the initial organic substrate concentration, the maximum voltage output decreased for low initial concentrations (0.125 g COD/L) and peaked for higher concentrations (0.25–4 g COD/L) at the same value (29 mV). Increasing the initial substrate concentration increased the duration of the maximum voltage. It was determined that an appropriate initial concentration for both high COD removal (90%) and lasting maximum voltage was approximately 2.5 g COD/L. Finally, when examining the effect of the electrolyte conductivity a maximum voltage was achieved (29 mV) by increasing the value beyond 2 S/m. These findings must be validated with respective experiments conducted using an H-type MFC. This model, despite the adequate MFC simulation, focuses on a specific MFC configuration with various assumptions. Further examination is needed in order to optimize and enrich the model with more electrochemical phenomena and obtain a more complete MFC simulation.