The Performance of a Microbial Fuel Cell Depends ... - Semantic Scholar

Bull Math Biol (2012) 74:834–857 DOI 10.1007/s11538-011-9690-0 O R I G I N A L A RT I C L E

The Performance of a Microbial Fuel Cell Depends Strongly on Anode Geometry: A Multidimensional Modeling Study Brian V. Merkey · David L. Chopp

Received: 10 March 2011 / Accepted: 5 August 2011 / Published online: 20 October 2011 © Society for Mathematical Biology 2011

Abstract A multidimensional biofilm model is developed to simulate biofilm growth on the anode of a Microbial Fuel Cell (MFC). The biofilm is treated as a conductive material, and electrons produced during microbial growth are assumed to be transferred to the anode through a conductive biofilm matrix. Growth of Geobacter sulfurreducens is simulated using the Nernst–Monod kinetic model that was previously developed and later validated in experiments. By implementing a conduction-based biofilm model in two dimensions, we are able to explore the impact of anode density and arrangement on current production in a MFC. Keywords Biofilm model · Microbial fuel cell 1 Introduction A Microbial Fuel Cell (MFC) is a promising technology that harnesses the natural growth mechanisms of particular bacterial species in order to produce electrical current. Bacterial growth (like growth of organisms in general) relies on the transfer of electrons from an initial electron donor to an eventual electron acceptor in order to produce energy (Rittmann and McCarty 2001). Bacteria growing aerobically use oxygen as an electron acceptor, while bacteria growing anaerobically will use a variety B.V. Merkey () · D.L. Chopp Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL, USA e-mail: [email protected] D.L. Chopp e-mail: [email protected] B.V. Merkey Mathematics Department, Viterbo University, La Crosse, WI, USA e-mail: [email protected]

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

835

of electron acceptors, including some solid materials. A MFC exploits those bacteria that naturally oxidize solid surfaces as part of their growth process, with the electrons deposited by these bacteria producing an electric current. Bacteria so exploited are sometimes referred to as anode-respiring bacteria (Rittmann et al. 2008). Though MFCs are not a recent discovery (Berk and Canfield 1964), there has in the past decade been a surge in research into MFCs (Logan et al. 2006; Logan and Regan 2006a). The reason for these efforts are many, from the renewable nature of energy produced from MFCs (Logan and Regan 2006b), to their use in low-power environmental sensors (Tender et al. 2002), as a natural source of hydrogen (Logan et al. 2008), and for their ability to convert waste material into usable energy (He et al. 2005; Min et al. 2005; Clauwaert et al. 2008; Lovley 2008; Rittmann et al. 2008). For all these reasons, the effort going into improving MFC systems and processes is certainly warranted; however, the success or failure of a MFC depends on many factors, many of which are poorly understood. In order to better understand MFC behavior, researchers have begun to turn to mathematical modeling as a tool to gain insight into these complex systems. Mathematical modeling of biofilms has been proven to be a successful technique for improving wastewater treatment processes (Rittmann and McCarty 2001; Wanner et al. 2006), and modeling of MFC processes has built on this success (Zhang and Halme 1995; Marcus et al. 2007; Picioreanu et al. 2007). Such modeling has, for example, shed light on the likely mechanisms employed by bacteria in transferring electrons to a MFC anode: for some organisms, the only plausible means of electron transfer is via direct conduction of electrons rather than via some diffusible electron shuttle or mediator (Torres et al. 2010). However, multidimensional modeling of the impact of electrode geometry and other operating conditions on the power output of a MFC has only been carried out for bacteria transferring electrons via a diffusible mediator (Picioreanu et al. 2007, 2010). Because there is evidence suggesting that direct conduction of electrons to a MFC anode is a more likely mechanism of electron transfer for an organism like Geobacter sulfurreducens (Reguera et al. 2006; Torres et al. 2010), in this work we develop a multidimensional model of biofilm growth in a MFC, where the biofilm is treated as a conductor. We base our model on the first conduction-based model of biofilm growth in a MFC (Marcus et al. 2007). This model has shown success in both explaining electron transfer in biofilms (Torres et al. 2008a, 2008b, 2010) and in measuring kinetic parameters for bacterial growth in a MFC (Lee et al. 2009). However, the 1-dimensional nature of this model carries with it the limitations of all 1D biofilm models (Wanner et al. 2006), and precludes easy study of the impact of electrode geometry or fluid flow conditions on biofilm development. This 1D model has since been improved to more accurately model recent experimental findings (Torres et al. 2008b; Marcus et al. 2010, 2011), but as a first step and to simplify model analysis we limit our model to the major features of the original 1D conduction-based model. This model is similar to previous multidimensional continuum models for biofilm growth (Alpkvist and Klapper 2007a; Duddu et al. 2008; Merkey et al. 2009). With our new model, we seek to explore the dependence of MFC power output on electrode geometry and other operational parameters when electron transfer occurs through a conductive biofilm. Note, though, that this model is intended to elucidate

836 Table 1 Unknown quantities solved for in the model

B.V. Merkey, D.L. Chopp Symbol Description

Unit

η

Local electric potential

V

Si

Concentration of substrate i

mg COD/cm3

Xj

Volume fraction of biomass species j

–

u

Fluid velocity outside the biofilm

mm/s

Φf

Fluid velocity potential

mm2 /s

ub

Biomass advective velocity inside the biofilm cm/d

Φ

Biomass velocity potential

cm2 /d

the trends one would expect to see in a physical MFC system, and is not meant to make precise predictions regarding the power output of a particular reactor configuration. In addition, this model sets the stage for further studies of questions of interest in microbial ecology, such as the evolution of anode-respiring bacteria (ARB) or the competition between ARB and other niche species.

2 Model Description This multidimensional model for MFC growth builds heavily on the 1D MFC model previously introduced by Marcus et al. (2007). Following their work, we model only the anode and assume sufficient proton transport to the cathode to facilitate electron transport into the anode by the biofilm. We treat the biofilm as a conductor with conductivity κbio , and we assume that pH and electroneutrality are both maintained within the biofilm. It has been shown that this latter assumption does not apply to all MFCs (see, e.g., Franks et al. 2009), and the build-up of protons in the biofilm as electrons are removed will lead to a decrease in the local pH that will inhibit bacterial growth (Torres et al. 2008b). This issue is being explored by other research groups (Marcus et al. 2011), and in this work we focus on questions related to the geometry of biofilm growth in a MFC. Knowing that biofilm growth in real MFCs is limited by proton transport tells us that this model may over-predict current when proton transport is a limitation, such as in very thick biofilms where there is strong diffusion resistance in general. The trends of current production as a function of system parameters, however, are likely to remain unchanged, and so long as we are aware of the limitations of these assumptions we are able to draw general conclusions from our model regarding the growth of biofilms in MFCs. Table 1 lists the unknown quantities that are solved for in this model, and the techniques used to obtain these solutions are explained in the sections that follow. 2.1 Computational Domain There are three primary regions that make up the computational domain: liquid, biofilm, and solid anode (Fig. 1); these are labeled respectively as Ωliquid , Ωbiofilm , and Ωanode . Like in previous biofilm models (Alpkvist and Klapper 2007b; Duddu et al. 2008; Smith et al. 2007; Merkey et al. 2009; Vaughan et al. 2010), the

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

837

Fig. 1 The computational domain with a summary of the equations being solved in each region

irregularly-shaped interfaces defining these regions (Γsurface and Γanode ) are described using the Level Set Method (Osher and Sethian 1988). Briefly, the interface is tracked using the zero-level contour of some function ϕ( x ) that is updated in time following the processes that affect the interface shape (such as biomass growth and interface erosion due to external fluid shear). 2.2 Electric Potential Following Marcus et al. (2007), an electron balance between the biofilm and anode yields the following equation for steady-state current density: ∇ · j −



ξk · Rk = 0,

(1)

k

where j is the current density (mA/cm2 ), Rk is the kth kinetic reaction involving movement of electrons to or from the anode, and ξk is the yield coefficient of electrons for reaction Rk . This equation holds everywhere inside the biofilm region Ωbiofilm . Using Ohm’s law, we may relate the current density j to the electrical potential η (volt): κbio ∇η + j = 0,

(2)

where κbio is the biofilm conductivity (mS/cm), and η is the local potential (volt) with respect to a reference half-max-rate potential: η ≡ E − EKA . The value of EKA depends strongly on the composition of the biofilm being studied (Lee et al. 2009). The form of (2) is similar to Fick’s law for diffusion, and like diffusion we assume that conduction has no inherent directionality. Combining (1) and (2) yields an elliptic

838

B.V. Merkey, D.L. Chopp

equation for the potential η that will be solved everywhere in the biofilm: κbio ∇ 2 η +



ξk · Rk = 0.

(3)

k

At the anode surface Γanode , (3) has a constant boundary condition, η|Γanode = Vanode ,

(4)

under the assumption that the anode potential is maintained fixed with a potentiostat. Note that the constant Vanode , like the potential η itself, is set with respect to the half-max-rate anode potential EKA : Vanode ≡ E|Γanode − EKA . At the biofilm/liquid interface Γsurface we have a no-flux condition for η:  ∂η  = 0, (5) ∂n Γsurface which is imposed under the assumption that electrons are not conducted outside the biofilm by the bacteria (the anode is the only outlet for electrons produced in the biofilm). Using (2), this condition may also be written as j|Γsurface = 0. 2.3 Substrates The mass balance on substrate concentrations in the biofilm follows previous biofilm models (Wanner et al. 2006). Also, like previous models, we use the difference in process time scales to assume a steady-state profile for substrate concentrations in the domain on the time scale of interest: biomass growth (Picioreanu et al. 1999). Thus, for a particular biomass configuration the substrate profiles are in steady-state, but if the biomass configuration or flow field changes, the substrate profiles must be recalculated. In practice, this is done at every simulation time-step, as described in Sect. 2.9. Steady-state substrate concentrations are affected by diffusion, advection, and consumption or production by the biofilm, so for each substrate we must solve an equation of the form: u Si ) + DS,i ∇ 2 Si − ∇ · (



γi,k · Rk = 0,

(6)

k

where Si is the local concentration of substrate i, DS,i is the diffusivity of substrate i, u is the steady-state advective fluid velocity, γi,k is the yield coefficient of reaction Rk , and the sum is over all reactions affecting this solute. Note, though, that not all terms in (6) are included in all regions. In this model, we assume that all reactions are mediated by biomass, and so Rk is nonzero only inside the biofilm (Ωbiofilm ) because bacteria are present only there. We also assume that the liquid is essentially immobile inside the biofilm, and so u ≡ 0 for all points in Ωbiofilm . Also, in the biofilm, we multiply the default acetate diffusion coefficient by a factor of 0.2 to account for the fact that mass transfer by diffusion is lower in the biofilm than in open water

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

839

(Stewart 2003). At the biofilm/liquid interface Γsurface , we impose two conditions on (6) simultaneously: the first is a continuity of concentration condition [Si ]Γsurface = 0, and the second is a flux-balance condition   ∂Si DS,i = 0. ∂n Γsurface

(7)

(8)

([ · ] Denotes a difference taken across the interface.) The latter condition is imposed using the Immersed Interface Method (LeVeque and Li 1994). At the anode surfaces Γanode , we impose a no-flux condition (∂Si /∂n = 0) because the solid anode is impenetrable to liquid, and at the domain boundaries we assume a constant bulk-liquid concentration Si = SB , periodic conditions, or outflow conditions, depending on the particular configuration. The fluid velocity u is calculated assuming potential flow for the liquid (Pozrikidis 1997). This means that the fluid velocity is assumed to be the gradient of a potential: u = ∇Φf , and this potential Φf is found by solving Laplace’s equation everywhere in the fluid region along with corresponding inflow and outflow boundary conditions (Fig. 1). However, there may be cases when there is no open path for fluid motion from the inflow to the outflow boundary due to excess growth by the biofilm. When this occurs, u is set identically to zero, and instead a boundary layer is imposed outside the biofilm (Wäsche et al. 2002; Wanner et al. 2006): the substrate concentration is set to the bulk value (Si,B ) outside the boundary layer distance, and within the boundary layer (6) is solved using only the diffusive term. 2.4 Biomass Following Wanner and Gujer (1986) and many others since, biomass quantities are represented as volume fractions Xi , indicating the fraction of available space at a particular location that is occupied by species i; the volume fractions of all species sum to 1 at all locations and at all times:  Xi = 1. (9) biomass

This model includes two biomass types, active and inert, with active biomass representing all components produced during growth (including actively growing bacteria and extracellular polymers), and inert biomass representing any other non-active components (such as “dead” cells). We make an additional assumption that the composition of the active biomass type does not change, with continuous production and hydrolysis of extracellular polymers yielding a roughly constant ratio of growing to produced biomass types. A mass balance on biomass quantities yields the following equation describing the behavior of the ith volume fraction quantity:  ∂Xi + ∇ · ( ub Xi ) = βi,k · Rk , ∂t k

(10)

840

B.V. Merkey, D.L. Chopp

where ub is the biomass advective velocity, and βi,k and Rk are reaction coefficients and rates, respectively, just as in (6). Note the time derivative included explicitly in (10), which captures the change in biomass occurring over the timescale of interest (hours or days); all other processes are assumed to be steady-state with respect to biomass growth (Picioreanu et al. 1999). The advective biomass velocity ub is calculated by combining (9) and (10) to yield    (11) βi,k · Rk . ∇ · ub = biomass

k

If we then assume smooth potential flow for all biomass elements (Klapper et al. 2002), we may write the biomass velocity as the gradient of a velocity potential Φ: ub = ∇Φ, and then (11) becomes:     ∇ 2Φ = βi,k · Rk . (12) biomass

k

This elliptic equation is solved only inside the biofilm region Ωbiofilm , for that is the only region where biomass movement occurs. The boundary conditions for (12) are no-flux at the anode surface Γanode because biomass may not move into the anode, a fixed value Φ = 0 at the biofilm/liquid interface Γsurface , and periodic on domain boundaries where appropriate. There are no boundary conditions imposed on nonperiodic domain boundaries because the simulations are stopped if the biofilm/liquid interface Γsurface reaches such a boundary. After Φ has been computed, the advective biomass velocity ub is found from ub = ∇Φ. Note that this advective velocity does not affect substrate concentration profiles within the biofilm because biomass advection is assumed to be slow compared to substrate diffusion (again see the separationof-timescales discussion in Sect. 2.3). 2.5 Evolution of the Biofilm/Liquid Interface As the biofilm grows, the location of Γsurface must be updated to account for the increased presence of biomass. Because the interface moves due to biomass movement, the normal interface velocity uint is calculated using the biomass velocity potential Φ:  ∂Φ  uint = . (13) ∂n Γsurface This velocity is then extended to all points in the domain using the Fast Marching Method (Adalsteinsson and Sethian 1995; Chopp 2001). In addition to interface movement due to biomass growth, we also include erosion of the biofilm surface. Like motion due to biomass growth, erosion forces move the interface with an erosion speed Fdet that in principle may take any form (Xavier et al. 2005). However, the physically-meaningful forms of Fdet are more limited (Morgenroth and Wilderer 2000; Picioreanu et al. 2001; Chambless and Stewart 2007), and we follow Duddu et al. (2009) in modeling erosion using the shear stress on

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

841

Fig. 2 Calculation of uint when Γsurface and Γanode intersect

the biofilm  surface due to fluid motion. Given a flow field and fluid/solid interface Γsurface Γanode , we calculate the tangential shear stress imposed on that interface by the fluid as  

∂u ∂v ∂v ∂u σ = μ n2x − n2y + + 2μnx ny − , (14) ∂y ∂x ∂y ∂x where μ is the kinematic viscosity of the fluid, nx and ny are the components of the interface normal vector, n = ∇ϕ/||∇ϕ||, and u and v are the components of the fluid velocity vector u . From the shear stress σ , we calculate an interface erosion speed Fdet = kDet |σ/σ1 |b , where kDet is an erosion strength coefficient with units of length per time, σ1 is a unit shear stress in units force per area, and b is a unitless exponent. Like the velocity of interface growth uint , the erosion velocity Fdet is extended outward from the interface to the rest of the domain using the Fast Marching Method. In this work, we use a fixed, assumed value for kDet (Sect. 2.8 and Table 3) and fix b at 0.58 based on previous literature (Rittmann 1982; Duddu et al. 2009). 2.6 Additional Interface Treatment for Biofilms Growing onto Solid Surfaces When a biofilm colony grows such that Γsurface intersects or grows into Γanode , extra treatment is required to maintain a well-behaved level set function ϕ( x ). We operate under the assumption that a biofilm attached to a solid surface will remain attached to this surface indefinitely, unless the biofilm/liquid interface Γsurface retreats such that the biofilm pulls away. Therefore, if the biofilm is attached to the surface (Γsurface lies parallel to and on top of Γanode ), we no longer use (13) to calculate the velocity of the interface because any perturbations or inaccuracies in this calculation can lead to artificial movement of the biofilm interface away from the solid surface. Instead, we consider a multipart calculation of uint (Fig. 2) that allows movement of the biofilm/liquid interface Γsurface along the anode surface Γanode , while also ensuring that an interface won’t artificially pull away from the anode surface due to slight inaccuracies in calculation of the interface velocity.

842

B.V. Merkey, D.L. Chopp

Table 2 Stoichiometric matrix of reaction kinetics Process

Acetate utilization

Biomass

Acetate

Potential

Xa

SA

η

Xi

− Y1 ρa

1

Biomass inactivation

−1

Biomass respiration

−1

ρa ρi

a

Kinetic equation

S

− Y1 λS ρa

1 A μa K +S Xa A A 1+exp(−(F /RT )η) bi Xa

−λX ρa

1 bres 1+exp(−(F /RT )η) Xa

a

Points on Γsurface outside Γanode are always treated as usual (uint = ∂Φ/∂n). However, points on Γsurface that lie inside Γanode are treated differently depending on their distance inside Γanode . Points on Γsurface that lie further inside Γanode than a prescribed distance use a zero-velocity condition: uint = 0. For points on Γsurface that lie closer to Γanode than this prescribed distance, we use a linear interpolation of the typical interface velocity and a zero velocity: uint = α ·

∂Φ ∂Φ + (1 − α) · 0 = α · , ∂n ∂n

where α is the distance inside Γanode as a fraction of the prescribed distance. In our model, we use one grid element as the prescribed distance over which this interpolation occurs, and so α = 0 for points on Γanode and α = 1 for points that are at a distance one grid element into the anode. 2.7 Kinetic Equations The kinetic equations used in this model are taken from Marcus et al. (2007) in a slightly rewritten form (Table 2). This model was well supported in later experimental work (Torres et al. 2008b; Lee et al. 2009), though for purposes of calculating the current density we retain the Nernst–Monod kinetic from Marcus et al. (2007) that was shown to be unnecessary for modeling bacterial growth by Torres et al. (2008a). 2.8 Model Parameters Parameters used in model simulations are given in Table 3. The model organism serving as the basis of active bacteria parameter values is Geobacter sulfurreducens, given its strong capacity for current production in MFCs (Bond and Lovley 2003; Reguera et al. 2006; Torres et al. 2010). The key growth parameters for G. sulfurreducens were determined by fitting experimental data from Bond and Lovley (2003) (Sect. 3.1), and the values were verified by comparison with other reported values or values used in previous modeling work (Marcus et al. 2007; Picioreanu et al. 2007; Lee et al. 2009). In order to determine a nonextreme erosion coefficient kDet , we simulated 10 days of growth under varying erosion strengths for flow speed 100 µm/s (Fig. 3). We chose the value kDet = 1.75 µm/d from the middle range of model outputs in order to avoid extreme effects. This value is the same order of magnitude as the erosion strength reported in Rittmann (1982).

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

843

Table 3 Default parameter values used in the model Parameter Description

Value

Units

Reference

μa

Specific growth rate of active biomass

0.69

d−1

F

Ya

Active biomass growth yield

0.049

g VS/g COD

F

KA

Half-max-rate acetate concentration

0.003

mg COD/cm3

F

bi

Biomass inactivation rate

0.08

d−1

F

bres

Biomass respiration rate

0.08

d−1

F A A

ρa

Density of active biomass

50

mg VS/cm3

ρi

Density of inactive biomass

50

mg VS/cm3

λS

Substrate-to-current conversion

0.13

mA · d/mg COD 1

λX

Biomass-to-current conversion

0.198

mA · d/mg VS

1

EKA

Half-max-rate anode potential

−0.448

V

1 2

Vanode

Fixed anode potential (Vanode = E|Γanode − EKA ) 0.648

V

κbio

Biofilm conductivity

0.5

mS/cm

3

DA

Diffusivity of acetate in water

0.941

cm2 /d

4 A

kDet

Erosion strength coefficient

1.75

µm/d

U

Influent flow speed

0.1

mm/s

A

Lb

Boundary layer thickness

50

µm

5

SB

Bulk acetate concentration (0.3 mM)

0.0191

mg COD/cm3

1

F

Faraday constant

96485

C/e−

4

R

Ideal gas constant

8.3145

J/mol/K

4

T

Temperature

303.15

K

2

μ

Kinematic viscosity of water

1.1 × 10−3 kg/m/s

4

Sources for values are: 1: Marcus et al. (2007), 2: Bond and Lovley (2003), 3: Torres et al. (2008a), 4: Lide (1990), 5: Lee et al. (2009), F: data fitting (Bond and Lovley 2003), A: assumed

Fig. 3 Dependence of biofilm area on erosion coefficient kDet . Simulations were run for 10 days with flow speed 100 µm/s and with varying only the erosion coefficient kDet . Low erosion values yield higher biomass production that is growth-limited, while for high erosion values biofilm growth is limited by erosion effects. Our chosen value of kDet = 1.75 µm/d lies between these behavioral extremes

844

B.V. Merkey, D.L. Chopp

2.9 Solution Strategy The algorithm used to model the growth of the biofilm over time is as follows. Note that the global time-step of simulations is set by the kinetics of biomass growth, which occur on the order of hours or days in comparison to the faster dynamics of substrate advection and diffusion and electron transport (Picioreanu et al. 1999). 1. Compute the fluid flow u around the biofilm and anode regions (Ωbiofilm ∪ Ωanode ). 2. Solve (3) and (6) for the electric potential and substrate profiles simultaneously by discretizing each equation and iterating with Newton’s Method until the profiles have converged. The Jacobian matrix required in the solution is made up of the finite-difference coefficients and analytically-calculated partial derivatives of the analytic reaction terms. 3. Compute the biomass growth velocity ub using (12), and update the biomass volume fractions using a conservative upwinding discretization of (10). 4. Use the Fast Marching Method to extend the biomass values outside the biofilm/ liquid interface Γsurface (Adalsteinsson and Sethian 1995; Chopp 2001), and update the location of Γsurface using the Level Set evolution equation (Osher and Sethian 1988). 5. Do a second update of Γsurface to capture erosion effects, with the interface velocity computed using the local fluid shear strength (14). For reasons related to computational convenience, this is carried out separate from the previous step. 6. Remove any portions of the biofilm that are no longer connected to the anode surface Γanode . This is accomplished by labeling points in the biofilm that neighbor the anode as connected, labeling neighbors of these points as also connected, and so on; any biofilm points that are not so marked must not have a connection to the anode, and will thus be removed (Xavier et al. 2005; Merkey et al. 2009). 7. Return to step 1, or quit if the simulation end time has been reached.

3 Results and Discussion In this work, we first validate our model by comparing with experimental results reported by Bond and Lovley (2003) (Sect. 3.1). Following model validation, we explore the effects of anode size and placement (Sect. 3.2), flow speed (Sect. 3.3), and anode density (Sect. 3.4) on current production in a MFC. When evaluating the performance of our model MFC, we consider the magnitude of total current density produced by the MFC. We denote this magnitude of current density by J , with units mA/cm2 , and J may be found by integrating the local current density over the entire anode surface:

1 J= j · nˆ dA, Aanode Γanode where nˆ is the unit normal into Γanode . The equivalent integral for the biofilm/liquid interface Γsurface evaluates to zero from (5), and we may add it to the above expression

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

845

in order to invoke Green’s Theorem and transform to an integral over the biofilm region Ωbiofilm :  

1 J= j · nˆ dS + j · nˆ dS Aanode Γanode Γsurface

1 ∇ · j dA. = Aanode Ωbiofilm Finally, making use of (1), we have the following expression for the net current produced by the biofilm:  

1 ξk · Rk dV J= Aanode Ωbiofilm k∈η

 −1 μa SA = λS + λX bres Aanode Ωbiofilm Ya KA + SA  1 × ρa Xa dV . 1 + exp(−(F /RT )η)

(15)

Table 2 illustrates that the local potential η will decrease in value as electrons are produced during biofilm growth, and the resulting gradient pushes electrons into the anode. However, by convention the flow of electrons is actually opposite to the direction of current, and this is the reason for the negative sign in (15): current actually flows into the biofilm, and so the outward flow is negative. In our results, we present current as a positive quantity with the understanding that all electron flows are into the anode surface. 3.1 Model Validation: Matching Bond and Lovley (2003) To validate our model and set parameter values to be used in the studies that follow, we matched our model output with experimental results from Bond and Lovley (2003). We inoculated a semicircular biofilm colony with radius 3 µm in an x-periodic domain of width 30 µm, and in order to recreate erosion effects due to fluid motion we limited biofilm growth to 8.5 µm above the solid anode surface. The remaining parameters are given in Table 3. In their experiments, Bond and Lovley (2003) reported that the bulk acetate concentration was replenished after 4 days and again after 5 days, with no replenishment from that point onward. Figure 4 illustrates the strong fit with experimental results we were able to obtain with our model. Because our model considered only current production from acetate as an electron donor, current production drops to zero when acetate has been exhausted (e.g., current production is zero at 4 days and again at 7 days in Fig. 4). However, Bond and Lovley reported a basal current production in acetate-depleted cultures that they determined was due to slow oxidation of nutrients present in the inoculum (such as succinate to fumarate); this basal production may be seen in the nonzero experimental values illustrated in Fig. 4. We included only acetate as an

846

B.V. Merkey, D.L. Chopp

Fig. 4 Comparison of model output with experimental results. Experimental data points from Bond and Lovley (2003) are plotted with circles, and direct model output is given by the solid line. An estimate of the nonzero basal current production reported by Bond and Lovley is plotted in the dotted line, and the dashed line illustrates model output for acetate-based current production combined with the basal rate

electron donor in our simulations and, therefore, in order to validate with these experimental results, we estimated their observed basal rate and used this rate combined with model output in our match with experimental data. The combined sources of current production yield a strong fit with experimental data, suggesting that a conduction-based biofilm model is appropriate for modeling growth of G. sulfurreducens. In validating their mediator-based MFC model, Picioreanu et al. (2007) predicted that a conduction-based model would likely yield a better fit with experiment than they were able to achieve, and we have shown that this is indeed the case. Parameter values resulting from the validation are summarized in Table 3. Encouragingly, our fit was obtained using parameter values similar (and in some cases identical) to those obtained by Picioreanu et al. (2007). 3.2 Effect of Anode Numbers When Anode Surface Area is Constant Having validated the model, in our first study, we sought to determine the effect of anode structure on biofilm growth and resultant current output. We varied the number of anode structures supporting biofilm growth while keeping the total anode surface area (2D perimeter) constant. Shown in Fig. 5 is the growth of biofilms under the conditions given in Table 3, but on different numbers of anode surfaces. At the inflow, the concentration of acetate is fixed at its bulk value (Table 3), and the concentration varies in space due to advection, diffusion, and consumption by the biofilm. The biofilm is inoculated solely with active biomass, and the single-anode case is inoculated to a uniform thickness of 20 µm; this is meant to simulate inoculation of the MFC in batch mode (Bond and Lovley 2003; Torres et al. 2008a; Lee et al. 2009). The simulations with higher anode numbers were inoculated with a decreased thickness that kept the total biofilm area constant between configurations. In all cases illustrated in Fig. 5, the biofilm colonies on each anode grow large enough to join with colonies from other anodes, resulting in one large biofilm colony. The predominant growth direction of this larger colony is toward the nutrient source at the left boundary. The biofilm has grown rather thick by 30 days, on the order of

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

847

Fig. 5 Biofilm growth in time for 1, 3, and 6 anodes. The left plots show the change in the biofilm interface location over time, and the right contour plots show acetate concentration after 30 days of growth. There is an inflow condition on the left domain boundary, an outflow condition on the right domain boundary, and periodic conditions for the upper and lower boundaries. The black circles represent the anode surfaces, and in the contour plots the white line indicates the location of the biofilm/liquid interface Γsurface

200 µm thick on the leading anode surfaces. In the 3 and 6-anode simulations, the biofilm has not grown to encompass the interior of the anode configuration, which is seen most easily in the time plots of interface locations in Fig. 5: the contours show almost no interface motion inside the anode structure. In all cases, the periodic domain has allowed the biofilm to grow across the region, leading to the possibility of flow channels being plugged (as has already happened in the 3-anode simulation). While it is possible that this is an artifact of simulations in a 2D domain, clogging by biomass is always a concern in biofilm systems in general (Logan et al. 2007;

848

B.V. Merkey, D.L. Chopp

Fig. 6 Current density over time (top) and after 25 days (bottom) for different anode numbers

Picioreanu et al. 2010), and reactor design generally takes this into account (He et al. 2005). The biofilm thicknesses obtained in these simulations are an order of magnitude larger than the thickness reported by, e.g., Lee et al. (2009). However, Marcus et al. (2007) reported that a conductivity κbio ≥ 10−3 mS/cm was sufficient to allow microbial growth tens of micrometers away from the anode surface; of course, the later estimate of κbio ≥ 0.5 mS/cm (Torres et al. 2008a) is well above this value. We also note that the fluid flow rate in this example (0.1 mm/s) is two or more orders of magnitude lower than is sometimes used in MFC experiments (see, e.g., Bond and Lovley 2003; Torres et al. 2009), resulting in lower overall fluid shear forces than may be seen experimentally. For these reasons, the development of a thicker biofilm in this example is not surprising. The sample images in Fig. 5 illustrate three cases out of several we considered, but current densities for all cases considered are shown in Fig. 6. All multianode cases showed transient behavior as biofilm growth led to closing of channels by the upstream anodes, leading to decreased current flux in the downstream anodes and

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

849

Fig. 7 Comparison of narrow (upper left), wide (upper right), and cross (lower left) anode positioning, along with current density over time (lower right)

a slight drop in overall current density. Between 10 and 25 days all anode configurations behave similarly, with a linear increase in current over time. By 30 days, though, the 3-anode configuration shows a drop in current production. This drop is due to closing of the flow channel, which means only the leading edge of the biofilm will be involved in current production rather than the entire biofilm surface. The other cases behave similarly over a longer period of time. The lower plot in Fig. 6 illustrates that the dependence of current density on anode number is minimal, with a variation of only 7% between all cases. Thus, we conclude that the actual surface area available for current production is more important than the shape configuration of the surface. This effect is likely a result of the high conductivity of the biofilm. The anode positions used in the previous simulations were all spaced regularly on the edge of the single-anode circle, but we also considered other distributions of anodes (Fig. 7). The purpose of these simulations was to test whether anode placement affects current production. Figure 7 illustrates the acetate contours after 30 days for three cases: 9 anodes arranged on a circle narrower than the single-anode circle, 9 anodes arranged on a larger circle, and 9 anodes arranged in a cross configuration. In addition, the current density is plotted as a function of time for these cases along with the 9-anode case plotted in Fig. 6. Initially, the wider anode spacing arrangement yields higher current production due to decreased mass transfer resistance of acetate to all anode surfaces, but over time the closing of the upstream channels leads to a drop in current production. For longer periods of time, all configurations show a similar progression of current production.

850

B.V. Merkey, D.L. Chopp

These examples all illustrate that maintaining sufficient interanode spacing to allow nutrient transport to all anode surfaces is crucial for achieving high performance in a MFC. 3.3 Effect of Flow Speed on Current Production We used the 6-anode configuration illustrated in Fig. 5 (bottom) to explore the dependence of current density on fluid flow speed. Shown in Fig. 8 is the dependence of current density on fluid flow speed. As the biofilms grow, current density generally increases over time, but for low flow speeds the rise in current is minimal due to mass transfer limitations. For intermediate flow speeds, there is generally a rise and then a drop in current, followed by a second rise; this results from closing of the upstream channels by the biofilm, limiting nutrient transport to the downstream surfaces. The highest flow speed erodes all biomass by around 10 days, sending current density to zero. Plotting the current density after 15 days as a function of flow speed (Fig. 8, bottom) illustrates that there is an optimal range of flow conditions to Fig. 8 Current density over time (top) and after 15 days (bottom) for varying flow speed

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

851

Fig. 9 Acetate concentration contours on 2, 8, and 32 anode bristles. Left: 10 µm-thick biofilm. Right: 25 µm-thick biofilm. The white line indicates the location of the biofilm/liquid interface Γsurface . All anode surfaces are contained within a 0.0196 mm2 region

maximize power output. Slow flows will not transport sufficient nutrients to support biofilm growth, while fast flows will shear excess biomass off the anode surfaces; a moderate flow avoids both these problems. Therefore, it seems important to adjust reactor recirculation rates to determine the appropriate balance of nutrient transport and erosion effects in MFC systems. Due to the scant experimental evidence on the response of a biofilm to flow conditions and the lack of a good estimate for the erosion parameter kDet , these simulations are meant as a qualitative rather than quantitative guide to reactor design. 3.4 Effect of Anode Density on Current Production Logan et al. (2007) first suggested using carbon fiber brushes as anode surfaces in a MFC. The bristles in these brushes are reported to have an average diameter of

852

B.V. Merkey, D.L. Chopp

Fig. 10 Current density (left) and current (right) for varying bristle density and biofilm thickness. In both plots, the anodes are contained within a 0.0196 mm2 region. The current density is calculated using (15), and the total current output by a particular configuration is calculated assuming a bristle length of 1.25 cm (Logan et al. 2007). For each biofilm thickness, mean values and standard deviations from 5 replicate simulations are plotted

7.2 µm, which is smaller than the anode sizes we have considered so far. In this section we explore the effect of bundling of these bristles on power output. We considered an increasing number of anode surfaces within a fixed geometrical area of (140 µm)2 = 0.0196 mm2 (Fig. 9), and considered several different biofilm thicknesses (5, 10, 15, and 25 µm) for each of these anode configurations. For each bristle density, we placed the anode surfaces randomly inside the fixed region, and carried out simulations for 5 such random placements to obtain a mean measure of model behavior. Initially, we consider only the immediate current production from a particular biofilm configuration rather than the development of the biofilm community in time. Currents output by various configurations are plotted in Fig. 10. Current density is greatest for low anode density, but the actual current produced by these configurations (current density times total anode surface area) peaks at intermediate to high anode densities. Higher current results from there being a large surface to support biomass growth but not so large to result in nutrient depletion, and these trends are a result of the interplay between the current production and nutrient consumption behaviors of the biofilm. In addition, thicker biofilms have more biomass available for current production and generally yield higher currents, but as anode numbers increase the biofilm reaches a saturation point where all available space is occupied by biomass (compare the left and right plots of Fig. 9). In this situation, increasing the biomass thickness or the number of anodes has a lower impact, and the biofilm may be viewed as a single large colony around a group of anodes rather than separate colonies on separate anodes. The variation in anode density may act as a proxy for distance from the brush center for the brushes described in Logan et al. (2007). Bristle density decreases with distance from the brush axis, so the plots in Fig. 10 suggest that current production may be higher at the ends of the bristles rather than nearer to the brush axis (where bristle density is highest). If it were not for the increased resistance that a longer bristle imposes, these results would suggest that longer bristles would yield higher power output than shorter bristles.

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

853

Fig. 11 Biofilm interface locations (left) and acetate concentration contours (right) for 2, 8, and 32 anode bristles. As the biofilm colonies grow, they merge and clog the flow field, changing the growth conditions

The variation in biofilm thickness yields some insight into the expected behavior of these configurations as the biofilm grows. When bristle density is high, the barriers to nutrient transport are high and increased biomass only makes them worse, so current density and total current will grow little as the biofilm grows. When bristle density is low, though, there are fewer barriers to nutrient transport, and increased biomass amounts will yield increased current production. Thus we expect that a lower anode density will yield the best anode performance in terms of current density. But, of course, we can also examine directly the growth in time of these different anode configurations. We considered growth in time for 2, 8, and 32 anode bristles in the domain, and again considered the mean behavior across 5 random anode placements. Similar to the behavior seen in Sect. 3.2, Fig. 11 illustrates that it is likely that a growing biofilm will bridge the gaps between bristles and merge to form one large

854

B.V. Merkey, D.L. Chopp

Fig. 12 Closeup of vectors of current production for thinner (left) and thicker (right) bristle biofilms. The white line indicates the location of the biofilm/liquid interface Γsurface , and the white vectors indicate the magnitude of current production on the anode surfaces. The current vectors tend to be shorter in the regions of lower acetate concentration, which are more pronounced for thicker biofilm colonies (right)

Fig. 13 Current density over time for varying bristle density. Plotted are the mean (symbols) and standard deviation (dashed lines) of current over time for 5 replicates of simulations with varying bristle numbers. Lower bristle densities yield higher current densities because there are fewer nutrient transport limitations for fewer anode surfaces. This echoes the results from Fig. 10

colony. The interior of this larger colony faces increased mass transfer resistance and a resulting decrease in current delivered to interior anode surfaces, as illustrated in Fig. 12. As a result, the current densities for higher bristle densities (which will clog more easily) tend to be lower than for lower bristle densities (which will be less prone to clogging), as shown in Fig. 13. Thus, our conclusions from carrying out simulations in time match the predictions made when considering only the biofilm community thickness. One final point of note is that, while it was shown in Fig. 10 that thicker biofilms tend to yield higher current production, the increased mass transfer resistance and lower current density seen in Fig. 12 suggest that the power output of a colony of this size is not as high as if this biomass were spread among a higher number of anode surfaces with more space between them. The implication of these results is that

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

855

maintaining sufficient nutrient delivery to all biofilm surfaces is crucial to maintaining high current production in a MFC.

4 Conclusions In this work, we have defined a framework for multidimensional modeling of a conductive biofilm in a MFC, and have shown that it is possible to extend the previous conductive-biofilm MFC model of Marcus et al. (2007) into multiple spatial dimensions. With multidimensional simulations, we are able to explore different questions than the previous work was designed to answer. Our simulations have shown that, at the anode scale, the primary driver of MFC current production is nutrient delivery to the entire biofilm surface. If access to nutrients is limited for any reason, then current production will suffer. This was seen most notably when biofilm colonies grew large enough to merge the anode array and limit nutrient transport to interior regions of the biofilm. It should be pointed out that the purpose of this work is to illustrate the trends one would expect to see in experimental measurements, and that the actual numbers reported by this model are not meant to be exact predictions of MFC performance. Acknowledgements We want to thank Bruce Rittmann for helpful discussions regarding this work. This work was supported by the NSF (Grant DMS-0921015) and the Initiative for Sustainability and Energy at Northwestern University (ISEN).

References Adalsteinsson, D., & Sethian, J. (1995). A fast level set method for propagating interfaces. J. Comput. Phys., 118(2), 269–277. Alpkvist, E., & Klapper, I. (2007a). Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction. Water Sci. Technol., 55(8–9), 265–273. Alpkvist, E., & Klapper, I. (2007b). A multidimensional multispecies continuum model for heterogeneous biofilm development. Bull. Math. Biol., 69, 765–789. Berk, R., & Canfield, J. (1964). Bioelectrochemical energy conversion. Appl. Environ. Microbiol., 12(1), 10. Bond, D. R., & Lovley, D. R. (2003). Electricity production by Geobacter sulfurreducens attached to electrodes. Appl. Environ. Microbiol., 69(3), 1548–1555. Chambless, J. D., & Stewart, P. S. (2007). A three-dimensional computer model analysis of three hypothetical biofilm detachment mechanisms. Biotechnol. Bioeng., 97(6), 1573–1584. Chopp, D. (2001). Some improvements of the fast marching method. SIAM J. Sci. Comput., 23, 230. Clauwaert, P., Aelterman, P., Pham, T. H., De Schamphelaire, L., Carballa, M., Rabaey, K., & Verstraete, W. (2008). Minimizing losses in bio-electrochemical systems: The road to applications. Appl. Microbiol. Biotechnol., 79(6), 901–913. Duddu, R., Bordas, S., Chopp, D. L., & Moran, B. (2008). A combined extended finite element and level set method for biofilm growth. Int. J. Numer. Methods Eng., 74(5), 848–870. Duddu, R., Chopp, D. L., & Moran, B. (2009). A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment. Biotechnol. Bioeng., 103(1), 92–104. Franks, A., Nevin, K., Jia, H., Izallalen, M., Woodard, T., & Lovley, D. (2009). Novel strategy for threedimensional real-time imaging of microbial fuel cell communities: monitoring the inhibitory effects of proton accumulation within the anode biofilm. Energy Environ. Sci., 2(1), 113–119. He, Z., Minteer, S., & Angenent, L. (2005). Electricity generation from artificial wastewater using an upflow microbial fuel cell. Environ. Sci. Technol., 39(14), 5262–5267.

856

B.V. Merkey, D.L. Chopp

Klapper, I., Rupp, C. J., Cargo, R., Purvedorj, B., & Stoodley, P. (2002). Viscoelastic fluid description of bacterial biofilm material properties. Biotechnol. Bioeng., 80(3), 289–296. Lee, H., Torres, C., & Rittmann, B. (2009). Effects of substrate diffusion and anode potential on kinetic parameters for anode-respiring bacteria. Environ. Sci. Technol., 43(19), 7571–7577. LeVeque, R. J., & Li, Z. (1994). The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Numer. Anal., 31(4), 1019–1044. Lide, D. R. (1990). CRC handbook of chemistry and physics (70th edn.). Boca Raton: CRC Press. Logan, B., Call, D., Cheng, S., Hamelers, H., Sleutels, T., Jeremiasse, A., & Rozendal, R. (2008). Microbial electrolysis cells for high yield hydrogen gas production from organic matter. Environ. Sci. Technol., 42(23), 8630–8640. Logan, B., Cheng, S., Watson, V., & Estadt, G. (2007). Graphite fiber brush anodes for increased power production in air-cathode microbial fuel cells. Environ. Sci. Technol., 41(9), 3341–3346. Logan, B. E., Hamelers, B., Rozendal, R., Schroder, U., Keller, J., Freguia, S., Aelterman, P., Verstraete, W., & Rabaey, K. (2006). Microbial fuel cells: Methodology and technology. Environ. Sci. Technol., 40(17), 5181–5192. Logan, B. E., & Regan, J. M. (2006a). Electricity-producing bacterial communities in microbial fuel cells. Trends Microbiol., 14(12), 512–518. Logan, B. E., & Regan, J. M. (2006b). Microbial fuel cells-challenges and applications. Environ. Sci. Technol., 40(17), 5172–5180. Lovley, D. R. (2008). The microbe electric: Conversion of organic matter to electricity. Curr. Opin. Biotechnol., 19(6), 564–571. Marcus, A., Torres, C., & Rittmann, B. (2010). Evaluating the impacts of migration in the biofilm anode using the model PCBIOFILM. Electrochim. Acta, 55, 6964–6972. Marcus, A., Torres, C., & Rittmann, B. (2011). Analysis of a microbial electrochemical cell using the proton condition in biofilm (PCBIOFILM) model. Bioresour. Technol., 102(1), 253–262. Marcus, A. K., Torres, C. I., & Rittmann, B. E. (2007). Conduction-based modeling of the biofilm anode of a microbial fuel cell. Biotechnol. Bioeng., 98(6), 1171–1182. Merkey, B. V., Rittmann, B. E., & Chopp, D. L. (2009). Modeling how soluble microbial products (SMP) support heterotrophic bacteria in autotroph-based biofilms. J. Theor. Biol., 259(4), 670–683. Min, B., Kim, J., Oh, S., Regan, J., & Logan, B. (2005). Electricity generation from swine wastewater using microbial fuel cells. Water Res., 39(20), 4961–4968. Morgenroth, E., & Wilderer, P. A. (2000). Influence of detachment mechanisms on competition in biofilms. Water Res., 34(2), 417–426. Osher, S., & Sethian, J. (1988). Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys., 79(1), 12–49. Picioreanu, C., Kreft, J. U., Klausen, M., Haagensen, J. A. J., Tolker-Nielsen, T., & Molin, S. (2007). Microbial motility involvement in biofilm structure formation—a 3D modelling study. Water Sci. Technol., 55(8–9), 337–343. Picioreanu, C., van Loosdrecht, M., Curtis, T., & Scott, K. (2010). Model based evaluation of the effect of pH and electrode geometry on microbial fuel cell performance. Bioelectrochemistry, 78(1), 8–24. Picioreanu, C., van Loosdrecht, M. C., & Heijnen, J. J. (1999). Discrete-differential modelling of biofilm structure. Water Sci. Technol., 39(7), 115–122. Picioreanu, C., van Loosdrecht, M. C., & Heijnen, J. J. (2001). Two-dimensional model of biofilm detachment caused by internal stress from liquid flow. Biotechnol. Bioeng., 72(2), 205–218. Pozrikidis, C. (1997). Introduction to theoretical and computational fluid dynamics. Oxford: Oxford University Press. Reguera, G., Nevin, K., Nicoll, J., Covalla, S., Woodard, T., & Lovley, D. (2006). Biofilm and nanowire production leads to increased current in Geobacter sulfurreducens fuel cells. Appl. Environ. Microbiol., 72(11), 7345. Rittmann, B. E. (1982). The effect of shear stress on biofilm loss rate. Biotechnol. Bioeng., 24(2), 501–506. Rittmann, B. E., & McCarty, P. L. (2001). Environmental biotechnology: principles and applications. New York: McGraw-Hill. Rittmann, B. E., Torres, C. I., & Marcus, A. K. (2008). Understanding the distinguishing features of a microbial fuel cell as a biomass-based renewable energy technology. In V. Shah (Ed.), Emerging environmental technologies (pp. 1–28). Berlin: Springer. Smith, B. G., Vaughan, B. L., & Chopp, D. L. (2007). The eXtended Finite Element Method for boundary layer problems in biofilm growth. Commun. Appl. Math. Comput. Sci., 2(1), 35–56. Stewart, P. (2003). Diffusion in biofilms. J. Bacteriol., 185(5), 1485.

The Performance of a Microbial Fuel Cell Depends Strongly on Anode

857

Tender, L., Reimers, C., Stecher, H., Holmes, D., Bond, D., Lowy, D., Pilobello, K., Fertig, S., & Lovley, D. (2002). Harnessing microbially generated power on the seafloor. Nat. Biotechnol., 20(8), 821– 825. Torres, C., Krajmalnik-Brown, R., Parameswaran, P., Marcus, A., Wanger, G., Gorby, Y., & Rittmann, B. (2009). Selecting anode-respiring bacteria based on anode potential: Phylogenetic, electrochemical, and microscopic characterization. Environ. Sci. Technol., 43(24), 9519–9524. Torres, C., Marcus, A., Lee, H., Parameswaran, P., Krajmalnik-Brown, R., & Rittmann, B. (2010). A kinetic perspective on extracellular electron transfer by anode-respiring bacteria. FEMS Microbiol. Rev., 34(1), 3–17. Torres, C., Marcus, A., Parameswaran, P., & Rittmann, B. (2008a). Kinetic experiments for evaluating the Nernst–Monod model for anode-respiring bacteria (ARB) in a biofilm anode. Environ. Sci. Technol., 42(17), 6593–6597. Torres, C., Marcus, A., & Rittmann, B. (2008b). Proton transport inside the biofilm limits electrical current generation by anode-respiring bacteria. Biotechnol. Bioeng., 100(5), 872–881. Vaughan, B., Smith, B., & Chopp, D. (2010). The influence of fluid flow on modeling quorum sensing in bacterial biofilms. Bull. Math. Biol., 72, 1143–1165. Wanner, O., Eberl, H., van Loosdrecht, M., Morgenroth, E., Noguera, D., Picioreanu, C., & Rittmann, B. (2006). Mathematical modeling of biofilms. Technical report, International Water Association. Wanner, O., & Gujer, W. (1986). A multispecies biofilm model. Biotechnol. Bioeng., 28, 314–328. Wäsche, S., Horn, H., & Hempel, D. C. (2002). Influence of growth conditions on biofilm development and mass transfer at the bulk/biofilm interface. Water Res., 36(19), 4775–4784. Xavier, J., Picioreanu, C., & van Loosdrecht, M. C. (2005). A general description of detachment for multidimensional modelling of biofilms. Biotechnol. Bioeng., 91(6), 651–669. Zhang, X.-C., & Halme, A. (1995). Modelling of a microbial fuel cell process. Biotech. Lett., 17(8), 809– 814.