Modelling Microbial Fuel Cells using lattice Boltzmann methods
Abstract
An accurate modelling of bioelectrochemical processes that govern Microbial Fuel Cells (MFCs) and mapping their behaviour according to several parameters will enhance the development of MFC technology and enable their successful implementation in well defined applications. The geometry of the electrodes is among key parameters determining efficiency of MFCs due to the formation of a biofilm of anodophilic bacteria on the anode electrode, which is a decisive factor for the functionality of the device. We simulate the bioelectrochemical processes in an MFC while taking into account the geometry of the electrodes. Namely, lattice Boltzmann methods are used to simulate the fluid dynamics and the advection–diffusion phenomena in the anode compartment. The model is verified on voltage and current outputs of a single MFC derived from laboratory experiments under continuous flow.
keywords:
microbial fuel cells, lattice Boltzmann, modelling, agentbased model, bioelectrochemical processesIntroduction
Due to improvements in power and current outputs achieved in the last decades ieropoulos2010improved (), studies in Microbial Fuel Cells (MFCs) have intensified santoro2017microbial (). Despite the increased interest in this field, the majority of research work focuses on practical experimental analysis. On the other hand, just a handful of studies ortiz2015developments () have steered towards mathematically describing and simulating the bioelectrochemical processes that generate the electrical output of these devices. Arguably, the intrinsic complexity of these systems requires highly elaborative procedures to produce accurate models. However, the advantages (in time and cost saving terms) acquired by an accurate MFC model has motivated researchers korth2015framework (); recio2016combined (); tsompanas2017cellular () to exert more effort towards this direction.
The utilisation of porous electrodes in continuous flow MFCs, which is realised as best practice due to higher outputs kim2012porous (), and the observation that the macrostructure of the anode, its material and its inhomogeneities have a significant impact Merkey2012 (); gajda2015simultaneous (), prompted the employment of the lattice Boltzmann method (LBM) in the framework of the MFC model presented here. LBM is a numerical tool based on mesoscopic phenomena and is able to model fluid dynamics and transport (advection–diffusion) phenomena in fluid systems Chopard2002 (). In contrast to its precursor, lattice gas automata that studies the individual molecular dynamics of fluids or gases –in a microscopic scale–, LBM simulates transport phenomena by studying the evolution of ensembles of molecules –in a mesoscopic scale. Nonetheless, the basic principles of microscopic physics are retained due to the utilisation of the ensemble averaged distribution function as the primary variable Chopard2002 (). LBMs have an inherent parallel nature that is advantageous when implementing the algorithm in specialised parallel hardware. Moreover, the lattices used can trivially illustrate complex geometries compared with other computation methods.
The ultimate goal of the proposed model is to simulate in an efficient way the behaviour of an MFC device, namely the voltage and current output. That is achieved by analysing in the macroscopic and mesoscopic scales several phenomena occurring in the anode halfcell of an MFC –assuming that the cathode performance is not limiting. Specifically, fluid dynamics, transport phenomena, bioelectrochemical processes and biofilm formation are all taken into consideration and, thus, simulated by separate procedures, that constitute the proposed model whilst interacting with each other. The formation of a biofilm of anodophilic bacteria on the anode electrode is a decisive factor for the functionality of the device. Fluid dynamics and transport phenomena were simulated using two different LBMs, one providing the velocity of the fluid throughout the porous electrode positioned in the anode compartment and the other calculating the concentration of chemicals –given the velocity field acquired from the first LBM.
The use of LBMs, in addition to providing the ability to represent inhomogeneities or complex geometries of electrodes, empowers the faster execution of critical calculations that can be considered as the bottleneck of previous studies that have studied these phenomena picioreanu2010model (). That is due to the numerous iterations needed for the system to reach an equilibrium in the procedures of fluid dynamics and transport phenomena and the frequency at which these two procedures have to be updated, namely every time the biofilm formation is updated.
To support the aforementioned statement, note that the deciding chemical, physical and biological processes occurring in an MFC need to be studied across a wide span of time scales von2009three (); bottero2013biofilm (). Although the time for a system to reach an equilibrium concerning its fluid dynamics and transport of chemicals, is in the period of some seconds, the evolution –or even a small change of the formation– of a biofilm is noted after some hours or days. As a result, when modelling the evolution in the behaviour of a system that includes such a complicated collection of interdisciplinary processes, the faster relaxing (reaching equilibrium) procedures, here fluid dynamics and advectiondiffusion, will have to be revised every time the slower updating process provides slightly different outputs. That introduces a bottleneck in the efficient execution of such detailed models.
1 Previous work
LBMs have been extensively studied and claimed to be a viable alternative to other computational methods Chopard2002 (). This method can be used for different modelling purposes, such as reactiondiffusion systems, complex fluid dynamics or wave propagation Chopard2002 (); wolf2004lattice (). LBM is considering a mesoscale statistical representation of groups of molecules and is divided into two processes, the collision and the propagation of these groups. The collision operator characterising each LBM, leads to the relaxation of the distributions of the groups to a local equilibrium which is the result of a set of parameters like the velocity, the concentration (considered as conserved quantities) and external forces.
The robustness of different types and variations of LBMs have been studied for several problems. For instance, in pan2006evaluation (), the fluid dynamics through complex geometries, namely porous media, was evaluated for several types of LBM and with different boundary conditions. In prasianakis2013simulation (), a threedimensional LBM was suggested, implementing twenty seven characteristic velocities to model the flow field in porous media. The authors paralleled the results from their model with the ones from other numerical models and measurements from experiments.
In addition to fluid dynamics, LBM has successfully tackled the solution of advection–diffusion equations with high accuracy chopard2009lattice (). Moreover, the ability of LBM to model nonlinear convection–diffusion phenomena with anisotropic diffusion have been presented in shi2011lattice (), where the modelling of several one and twodimensional paradigms was investigated in detail. In li2017lattice () a comparison of different lattices has been executed to analyse the robustness and accuracy of them. Also, the impact that the chosen scheme of boundary conditions has in the accuracy of the model was studied huang2015boundary ().
Moreover, researchers have previously employed the LBM in the field of conventional fuel cells. For instance, in joshi2007lattice () the diffusion of gases –that are characterised by different molecular weights– were modelled in complex porous media, specifically in a solid oxide fuel cell (SOFC) anode. The authors in wang2006modeling () studied the flow field in a threedimensional scheme of a section of a serpentine channel and in a twodimensional scheme of a channel filled or partially filled with a porous medium, taking particular interest in the boundary conditions and the simulated time periods. Moreover, in dang2016pore () a multicomponent flow model was presented to analyse the transport of a mixture of gases in a SOFC anode through irregular porous media. The impact of the microstructure of the electrode, the dimensionless reaction flux and fuel composition to the fuel cell outputs were studied. Also, the effect that the microstructure of the anode has on the distribution and deposition of carbon and, thus, to the performance of the cell was studied xu2016lattice ().
In addition, LBMs have been used in threedimensional simulation of biofilm growth and evolution with respect to the flow dynamics and the advection–diffusion of chemicals through porous media von2009three (); pintelon2012effect (). Two types of biofilms were studied, with zero permeability von2009three () and nonzero permeability pintelon2012effect (). The main difference in the work presented here with the aforementioned, is that here the function of the anode compartment of an MFC is simulated, namely, the bioelectrochemical reactions were studied in addition to the flow dynamics, transport phenomena and biofilm formation. Nonetheless, the impact that the bioelectrochemical phenomena have towards all the other procedures were included to the model. Keeping in mind the ultimate goal of an enhanced performance of the model, the presented model is based on a unified discrete grid where all procedures of the model are calculated. On the other hand, previous works von2009three (); pintelon2012effect () have used LBMs based on discrete grids, whereas the simulation of biofilm is tackled with an individualbased model picioreanu2004particle () utilising continuous Cartesian dimensions. This incoherency in the procedures of the model will pose difficulties while implementing the model efficiently in software (or indeed hardware).
In the present work the selection of using a twodimensional grid, rather than a threedimensional similar to previous studies, certainly produces a lower accuracy, as claimed in von2009three (); pintelon2012effect (). However, the use of a twodimensional grid enables a faster execution of the model, whilst scaling up the grid to three dimensions is trivial. Moreover, the utilisation of the same grid for all the procedures of the model, is an additional advantage when advancing towards three dimensions.
As mentioned before, a great advantage of LBMs is their capacity of significant acceleration which is straightforward, due to the inherently parallel characteristics of LBM. The parallel implementation of the algorithm of the model can be realised with advanced programming languages compatible with multicore processors bacsaugaouglu2016computational () or specialised hardware –like Field Programmable Gate Arrays (FPGAs) 4439254 (); 903392 () or Graphics Processing Units (GPUs) 5362489 ().
The biofilm formation is a process that it is not yet fully decoded by scientific research. The accurate modelling of biofilms is a quite demanding task, given that the behaviour of living, selforganising organisms has to be taken into consideration and the experimental justification of simulated results is challenging. However, several examples of previous work aimed to explain and model the underlying mechanics of the evolution of the biofilm, like attachment, growth, decay, lysis and detachment of biomass bottero2013biofilm (); picioreanu2004particle (); stewart2002pore (); bol20093d (). A simple agentbased model that can simulate the basic procedures of attachment and development was used here. Note here, that despite the simplicity of the model, it is capable of robustly formulating the under study processes.
Further advantages of the LBM model are that a detailed analysis of mass transfer can be carried out for the actual microstructure of the anode electrode and the cell’s performance can be judged by calculating the concentration polarisation. Complex geometries can be easily handled and both continuum and noncontinuum diffusion through the pores can be modelled. Because of the level of detail, the effect of anode geometry on voltage polarisation can be directly studied and the geometry can be optimised to provide enhanced performance. In fact, being able to model an anodic microstructure may be advantageous in understanding 3D porous electrode behaviour of real MFCs and their biofilms.
2 Methods
To approximate the behaviour of a MFC the model simulates the formation of a biofilm on the electrode surface, the fluid dynamics and the advection–diffusion phenomena throughout the anode compartment and calculates outputs of the MFC based on local concentration of chemicals and electrochemical equations. Two LBMs were used. One LBM represents the dynamics of fluids inside the anode compartment and another LBM represents advectiondiffusion equation of chemicals. The model runs as follows:

Formation of the electrode geometry and initialisation of the parameters.

Calculation of fluid dynamics with LBM (first LBM model).

Calculation of the advection–diffusion equation with LBM (second LBM model).

Calculation of MFC outputs (bioelectrochemical model).

Biofilm formation.

Go to step (2).
2.1 Lattice Boltzmann Method
A LBM is characterised by three main factors: the lattice, the equilibrium equation and the kinetic equation. We adopt the BhatnagarGrossKrook (BGK) LBM framework Qian92 () with a D2Q9 lattice (Fig. 1), based on a twodimensional lattice with nine characteristic velocities:
(1) 
(2) 
(3) 
The particle velocities distribution is governed by the evolution equation which takes into account propagation and collision processes:
(4) 
where is the dimensionless relaxation time, is a simulation time step, is a site in the lattice, and is an index of the characteristic velocity. is the equilibrium distribution function:
(5) 
where is a model constant relating to the lattice speed (), are defined in Eqs. (13), is the velocity flow, is the mass density and is a weight factor equal to:
(6) 
(7) 
(8) 
2.1.1 Fluid dynamics
Fluid dynamics are calculated by LBM with boundaries defined by the compartment physical boundaries (i.e. plastic container), the anode electrode geometry and the dynamic biofilm formation. The method is based on the BGK model, see Eqs. (1 – 8). This method is utilising bounceback boundary conditions Chopard2002 () and has a second order accuracy for fluid dynamics. The resulting mass density () and momentum density () determine the flow field as follows:
(9) 
(10) 
On the boundaries the bounceback principle dictates that the evolution equation is not calculated as in Eq. (4). However, the collision process is updated depending on the configuration of the boundaries and the available space. For instance, in the case of a boundary placed in the south direction of the lattice (Fig. 2) the following equations apply:
(11) 
2.1.2 Advection–diffusion equations
The advection–diffusion equations are solved using another BGK LBM, where the flow field and geometry of electrodes are used as inputs. The boundary conditions are considered as Dirichlet boundaries for the input and output of the domain (with a given concentration amount) and as Neumann for the solid boundaries (with a set zeroflux of concentration towards them) zhang2012 (). The evolution equation for the distribution function is similar to Eq. (4) (note the different notation of the distribution function used – rather than – to distinguish the two separate LBMs):
(12) 
where is the equilibrium distribution function, which is calculated as in Eq. (5) and are the weights defined by Eqs. (6)(8). Note that the velocity () is considered as a constant for the advection–diffusion simulation and provided by the first LBM procedure.
The resulting local concentration () which represents the concentration field is calculated in the following :
(13) 
2.2 Bioelectrochemical model
The mechanism for the transfer of electrons to the anode is assumed to incorporate the existence of an intracellular mediator and the interconnection of the biomass with the electrode by nanowires or by direct contact recio2016combined (); pinto2010two (). In the simulated experiments the total concentration of the intracellular mediator is regarded as constant and analogous to the biomass concentration. After calculating the equilibrium of fluid dynamics and the concentration of the chemical species, the output of MFC is studied using the following equations:
(14) 
(15) 
where is the total amount of intracellular mediator, is the reduced form, while is the oxidised form. is the mediator yield granted by the utilisation of the substrate, is the mediator molar mass, is the fraction of electrons provided during reduction of the intracellular mediator. is the Faraday constant, is the anode compartment volume, is the concentration of biomass. is the total current produced and, thus, flowing through the device and is the consumption rate of the substrate.
The double Monod equation is used to determine the consumption rate of the substrate:
(16) 
where is the concentration of the substrate, is the half saturation constant for substrate, is the half saturation constant for oxidised mediator and is the maximum consumption rate.
The MFC as an electrochemical system suffers from losses, defined elsewhere as polarisation or overpotential, caused by the rates of the reactions (activation overpotential), the ohmic resistances of the material used to build the MFC (ohmic overpotential) and the limitations posed on the mass transfer or kinetic phenomena (concentration overpotential). As a result, considering Kirchhoff’s voltage law and Ohm’s law the following equation correlates the overpotentials, the values of resistance and the produced current.
(17) 
where and are the external (load) and the internal resistances, is the maximum produced voltage (or the open circuit voltage), and are the concentration and activation overpotentials, respectively. The concentration overpotential based on the Nerst equation is defined by:
(18) 
where is the universal gas constant and is the temperature; while, the activation overpotential can be calculated by an approximation of the ButlerVolmer equation:
(19) 
where is the anode electrode surface area and is the exchange current density for mediator oxidation in reference conditions. Consequently, the output current can be estimated by:
(20) 
where is a constant which bounds the output current at low values of . While, Ohm’s law provides the voltage produced:
(21) 
2.3 Biofilm formation
The biofilm formation is an agent–based procedure simulating the attachment and growth of biomass. The attachment of biomass is simulated by a random selection in each time step of a predefined number of lattice cells () to change their biomass concentration to a predefined initial value (). The random selection of these cells is amongst the ones representing the front between the fluid domain and the anode electrode domain. Note that the initial value of the intracellular mediator in the cells selected for biomass attachment is set and equivalent to the biomass predefined initial value.
The growth of the biomass is following the double Monod equation providing the consumption rate (Eq. (16)). Along with the growth of biomass, an analogous increase of the intracellular mediator is assumed. The spreading of the biomass is initiated when the biomass concentration in a lattice cell exceeds a predefined threshold (). The direction towards the spreading of the biomass happens, is randomly selected. If the under study cell is surrounded by cells representing unavailable space (solid boundaries and existing amounts of biomass) a random walk is executed to find a free cell towards where a displacement of biomass can occur. When a lattice cell representing an available choice is found, each cell with an existing value of biomass, transports its biomass value to the neighbour indicated by the random walk. Finally, when an adjacent cell of the initial cell exceeding the predefined threshold is freed up, a fraction of the biomass concentration of the initial cell () is subtracted from its biomass concentration and it is transported to the currently free cell. In addition to the transport of biomass being displaced by the random walk, the amounts of intracellular mediators follow the same procedure.
Following the completion of the biofilm formation procedure, all the algorithmic steps of the model are repeated, starting with the calculation of fluid dynamics with LBM as formerly defined. The fluid flow is recalculated taking into account the updated biofilm as impermeable. Consequently, the advection and diffusion of the chemicals have to also be updated. Then, the electrical outputs are revised and, finally, the biofilm with the new concentrations and flow field is reformed. The model was implemented with Matlab R2016b.
3 Results
The applicability of the model can be verified with the comparison of the simulated results with measurements obtained in the laboratory. More specifically, here, the measurements obtained by an experiment with an innovative configuration of a single MFC cell, located onto a MFC system based on a building brick, were used. The MFC is comprised by the anode compartment accommodated in one cavity of the brick, the cathode compartment in another cavity and the ceramic semipermeable separator between them.
The influnet of the anode compartment is considered for the simulation and depicted in the following figures to enter from the left side of the axis. Consequently, the effluent is considered to exit the anode compartment from the right side of the axis of the illustrated model configuration.
The values of the parameters used in the aforementioned equations that describe the functionality of the MFC are given in Table 1. These values are extracted by measurements in the laboratory experiment or adopted by previous attempts of describing MFCs and biofilms with mathematical equations in the literature bottero2013biofilm (); pinto2010two ().
Parameter  Description  Value 

Anode compartment dimensions  
Anode electrode dimensions  
Model lattice dimensions  
Electrode porosity  0.874  
Fluid input velocity  
Kinematic viscosity  
Dimensionless relaxation time (LBM)  0.6706  
Diffusion coefficient  
Dimensionless relaxation time (LBM)  0.5036  
Input substrate concentration  410  
Total amount of intracellular mediator  0.05  
Intracellular mediator yield  0.5687  
Mediator molar mass  663400  
Electrons provided during reduction  2  
Maximum consumption rate  8.48  
Half saturation constant for substrate  20  
Half saturation constant for oxidised mediator  
Exchange current density  
Open circuit voltage  
External ohmic resistance (load)  360  
Internal ohmic resistance  360  
Predefined number of cells for attachment  200  
Initial biomass concentration  450  
Threshold biomass concentration  512.5  
Fraction of the biomass spreading  40%  

The procedures of calculating fluid dynamics and the advection–diffusion equation are running for enough time stems until an equilibrium is reached, while the bioelectrochemical model and the procedure of biofilm formation are updated once during each simulation interval. One iteration of the model corresponds to one hour of real time.
In Fig. 3 the simulated anode compartment configuration is depicted. The black lines on top and on the bottom of the figure represents the impermeable boundaries (chassis) of the anode compartment. The influent is entering the compartment from and the effluent is exiting from . Whereas, the porous medium of the anode electrode is illustrated with the black shapes that are repeated throughout the figure. Note that a random geometry of the anode electrode was adopted. Nonetheless, this selection imposes no limitation on the functionality of the model. The formation of the geometry of the anode electrode can be disorganised, lacking a regular geometry or order, or based on mathematical selforganising mathematical tools bandman2011using ().
The flow and concentration field in the anode compartment can be calculated without any biofilm formed on the anode electrode. The results of such calculations are illustrated in Figs. 4 and 5.
When a biofilm is developed on the anode electrode (illustrated in Fig. 6 as the gray spots attached on the black shapes that represent the porous electrode) then the flow and concentration fields are recalculated taking into consideration the new configuration. The results of the flow field LBM calculation are depicted in Fig. 7 and the results of the advectiondiffusion LBM are depicted in Fig. 8.
Finally, the simulation of the functionality of the under study MFC is carried out for a simulation of 72 hours and the results are presented in the following figures. In Fig. 9 the output voltage of the MFC is illustrated. The output voltage is asymptotically increasing during the first 60 hours and then reaches a steady state of approximately 330 mV that is in agreement with the laboratory based measurements. Moreover, the output current depicted in Fig. 10 is following the increase of the voltage curve and it reaches a steady value of 0.91 mA (similar with the laboratory experiment).
In Fig. 11 the ratio of the reduced and the oxidised mediator to the total amount of intracellular mediator are depicted through time. The concentration of reduced mediator (solid line) is following the evolution of the output voltage, whereas the concentration of the oxidised mediator (dashed line) is following a complete opposite asymptotic curve of the reduced mediator, as derived by Eq. (14).
4 Conclusions
We proposed a model of MFC that incorporate geometry of the electrodes. The model implements two LBMs to calculate the fluid flow field and the advection–diffusion phenomena of chemicals inside the anode compartment. Also, an agent–based model calculates the attachment and spreading of biomass on the electrode. The output voltage is calculated using bioelectrochemical processes.
An advantage of the model is that it is designed to consider the flow dynamics of the influent and the advection–diffusion phenomena of chemicals in comparably short time intervals (seconds and minutes). Thus, the model can be utilised in predicting transitions of the behaviour of a MFC from one state to another in this time span when the relevant parameter, namely influent velocity or fuel concentration is altered.
We adopted a random geometry of the anode electrode. This imposes no limitations on the model’s functionality. The formation of the geometry of the anode electrode can be disorganised, lacking a regular geometry or order, or based on mathematical selforganising mathematical tools bandman2011using (). Studying particular geometries will be a topic of future research. This could include analysis of the effect of specified geometries and microstructures on the development of biofilms in a MFC and its electrical and chemical outputs, e.g. geometries based data from Xray tomographic microscopy prasianakis2013simulation ().
The biofilm formed during modelled development of a MFC is assumed to impose a very high resistance to fluid dynamics, thus, the speed of fluid components throughout the biofilm was considered negligible. However, this is not an ideal representation of a real biofilm. More work would be required to take into account the permeability of biofilms and communication between cells in the matrix.
While, in the simulated experiments the total concentration of the intracellular mediator is regarded as constant and analogous to the biomass concentration, in reality, fractions of oxidised and reduced forms of the mediator can vary. This could be another aspect of future work.
Acknowledgment
This work was supported by the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 686585.
References
 (1) I. Ieropoulos, J. Greenman, C. Melhuish, Improved energy output levels from smallscale microbial fuel cells, Bioelectrochemistry 78 (1) (2010) 44–50.
 (2) C. Santoro, C. Arbizzani, B. Erable, I. Ieropoulos, Microbial fuel cells: From fundamentals to applications. a review, Journal of Power Sources 356 (2017) 225–244.
 (3) V. OrtizMartínez, M. SalarGarcía, A. De Los Ríos, F. HernándezFernández, J. Egea, L. Lozano, Developments in microbial fuel cell modeling, Chemical Engineering Journal 271 (2015) 50–60.
 (4) B. Korth, L. F. Rosa, F. Harnisch, C. Picioreanu, A framework for modeling electroactive microbial biofilms performing direct electron transfer, Bioelectrochemistry 106 (2015) 194–206.
 (5) D. RecioGarrido, M. Perrier, B. Tartakovsky, Combined bioelectrochemical–electrical model of a microbial fuel cell, Bioprocess and biosystems engineering 39 (2) (2016) 267–276.
 (6) M.A. Tsompanas, A. Adamatzky, I. Ieropoulos, N. Phillips, G. C. Sirakoulis, J. Greenman, Cellular nonlinear network model of microbial fuel cell, BioSystems 156 (2017) 53–62.
 (7) J. R. Kim, H. C. Boghani, N. Amini, K.F. AgueyZinsou, I. Michie, R. M. Dinsdale, A. J. Guwy, Z. X. Guo, G. C. Premier, Porous anodes with helical flow pathways in bioelectrochemical systems: The effects of fluid dynamics and operating regimes, Journal of Power Sources 213 (2012) 382–390.

(8)
B. V. Merkey, D. L. Chopp, The
performance of a microbial fuel cell depends strongly on anode geometry: A
multidimensional modeling study, Bulletin of Mathematical Biology 74 (4)
(2012) 834–857.
doi:10.1007/s1153801196900.
URL https://doi.org/10.1007/s1153801196900  (9) I. Gajda, J. Greenman, C. Melhuish, I. Ieropoulos, Simultaneous electricity generation and microbiallyassisted electrosynthesis in ceramic mfcs, Bioelectrochemistry 104 (2015) 58–64.
 (10) B. Chopard, A. Dupuis, A. Masselot, P. Luthi, Cellular automata and lattice boltzmann techniques: An approach to model and simulate complex systems, Advances in complex systems 5 (02n03) (2002) 103–246.
 (11) C. Picioreanu, M. C. van Loosdrecht, T. P. Curtis, K. Scott, Model based evaluation of the effect of ph and electrode geometry on microbial fuel cell performance, Bioelectrochemistry 78 (1) (2010) 8–24.
 (12) D. von der Schulenburg, T. Pintelon, C. Picioreanu, M. Van Loosdrecht, M. Johns, Threedimensional simulations of biofilm growth in porous media, AIChE Journal 55 (2) (2009) 494–504.
 (13) S. Bottero, T. Storck, T. J. Heimovaara, M. C. van Loosdrecht, M. V. Enzien, C. Picioreanu, Biofilm development and the dynamics of preferential flow paths in porous media, Biofouling 29 (9) (2013) 1069–1086.
 (14) D. A. WolfGladrow, Latticegas cellular automata and lattice Boltzmann models: an introduction, Springer, 2004.
 (15) C. Pan, L.S. Luo, C. T. Miller, An evaluation of lattice boltzmann schemes for porous medium flow simulation, Computers & fluids 35 (8) (2006) 898–909.
 (16) N. Prasianakis, T. Rosén, J. Kang, J. Eller, J. Mantzaras, F. Büichi, Simulation of 3d porous media flows with application to polymer electrolyte fuel cells, Communications in Computational Physics 13 (3) (2013) 851–866.
 (17) B. Chopard, J. Falcone, J. Latt, The lattice boltzmann advectiondiffusion model revisited, The European Physical Journal Special Topics 171 (1) (2009) 245–249.
 (18) B. Shi, Z. Guo, Lattice boltzmann simulation of some nonlinear convection–diffusion equations, Computers & Mathematics with Applications 61 (12) (2011) 3443–3452.
 (19) L. Li, R. Mei, J. F. Klausner, Lattice Boltzmann models for the convectiondiffusion equation: D2Q5 vs D2Q9, International Journal of Heat and Mass Transfer 108 (2017) 41–62.
 (20) J. Huang, W.A. Yong, Boundary conditions of the lattice boltzmann method for convection–diffusion equations, Journal of Computational Physics 300 (2015) 70–91.
 (21) A. S. Joshi, K. N. Grew, A. A. Peracchio, W. K. Chiu, Lattice boltzmann modeling of 2d gas transport in a solid oxide fuel cell anode, Journal of power sources 164 (2) (2007) 631–638.
 (22) L.P. Wang, B. Afsharpoya, Modeling fluid flow in fuel cells using the latticeboltzmann approach, Mathematics and Computers in Simulation 72 (2) (2006) 242–248.
 (23) Z. Dang, H. Xu, Pore scale investigation of gaseous mixture flow in porous anode of solid oxide fuel cell, Energy 107 (2016) 295–304.
 (24) H. Xu, Z. Dang, Lattice boltzmann modeling of carbon deposition in porous anode of a solid oxide fuel cell with internal reforming, Applied Energy 178 (2016) 294–307.
 (25) T. R. Pintelon, C. Picioreanu, M. van Loosdrecht, M. L. Johns, The effect of biofilm permeability on bioclogging of porous media, Biotechnology and bioengineering 109 (4) (2012) 1031–1042.
 (26) C. Picioreanu, J.U. Kreft, M. C. van Loosdrecht, Particlebased multidimensional multispecies biofilm model, Applied and environmental microbiology 70 (5) (2004) 3024–3040.
 (27) H. Başağaoğlu, J. Blount, J. Blount, B. Nelson, S. Succi, P. M. Westhart, J. R. Harwell, Computational performance of sequencel coding of the lattice boltzmann method for multiparticle flow simulations, Computer Physics Communications.
 (28) K. Sano, O. Pell, W. Luk, S. Yamamoto, Fpgabased streaming computation for lattice boltzmann method, in: 2007 International Conference on FieldProgrammable Technology, 2007, pp. 233–236. doi:10.1109/FPT.2007.4439254.
 (29) M. Gokhale, J. Stone, J. Arnold, M. Kalinowski, Streamoriented fpga computing in the streamsc high level language, in: Proceedings 2000 IEEE Symposium on FieldProgrammable Custom Computing Machines (Cat. No.PR00871), 2000, pp. 49–56. doi:10.1109/FPGA.2000.903392.
 (30) P. Bailey, J. Myre, S. D. C. Walsh, D. J. Lilja, M. O. Saar, Accelerating lattice boltzmann fluid flow simulations using graphics processors, in: 2009 International Conference on Parallel Processing, 2009, pp. 550–557. doi:10.1109/ICPP.2009.38.
 (31) T. L. Stewart, H. Scott Fogler, Porescale investigation of biomass plug development and propagation in porous media, Biotechnology and bioengineering 77 (5) (2002) 577–588.
 (32) M. Böl, R. B. Möhle, M. Haesner, T. R. Neu, H. Horn, R. Krull, 3d finite element model of biofilm detachment using real biofilm structures from clsm data, Biotechnology and Bioengineering 103 (1) (2009) 177–186.
 (33) Y. H. Qian and D. D´Humières and P. Lallemand, Lattice bgk models for navierstokes equation, EPL (Europhysics Letters) 17 (6) (1992) 479.

(34)
T. Zhang, B. Shi, Z. Guo, Z. Chai, J. Lu,
General bounceback
scheme for concentration boundary condition in the latticeboltzmann method,
Phys. Rev. E 85 (2012) 016701.
doi:10.1103/PhysRevE.85.016701.
URL http://link.aps.org/doi/10.1103/PhysRevE.85.016701  (35) R. Pinto, B. Srinivasan, M.F. Manuel, B. Tartakovsky, A twopopulation bioelectrochemical model of a microbial fuel cell, Bioresource technology 101 (14) (2010) 5256–5265.
 (36) O. Bandman, Using cellular automata for porous media simulation, The Journal of Supercomputing 57 (2) (2011) 121–131.