Pierre Belleville

Gérard Merlin

Julien Ramousse

Jonathan Deseure

Two-dimensional modelling of syntrophic glucose conversion in bioanodes for coulombic efficiency optimization

Keywords: Modelling, biomass segregation, bioanode, biofilm, Simulated Pareto front Nomenclature

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

D

Diffusion coefficient / m² s -1

D fact

Relative effective coefficient of diffusion / - ν

∆G
Local biomass advection speed in biofilm / m s -1 q max, gibbs Maximum Energy consumption rate / J mol_X j -1 s -1 q max, i Maximum consumption rate / mol_i mol_X j -1 s -1 r i Consumption rate / mol_i s -1 R Gas constant / J mol -1 K -1

Introduction

Biomass is considered as a renewable energy source of major interest for the production of fuel and chemicals [START_REF] Langeveld | Development Perspectives Of The Biobased Economy: A Review[END_REF]. From all, organics matter in wastewater have a great potential of valorization [START_REF] Rozendal | Towards practical implementation of bioelectrochemical wastewater treatment[END_REF]. Microbial ability to use a solid conductive material as the final electron acceptor to run their metabolism has led to the emergence of a quickly expanding research area, bioelectrochemical systems (BES) [START_REF] Bond | Electricity production by Geobacter sulfurreducens attached to electrodes[END_REF]. BES is a promising technology to couple wastewater treatment and valuable products (energy or chemicals) [START_REF] Logan | Conversion of wastes into biolectricity and chemicals using microbial electrochemical technologies[END_REF]. This principle can be applied to different BES applications such as power production [START_REF] Logan | Microbial Fuel Cells: Methodology and Technology[END_REF], hydrogen production [START_REF] Logan | Microbial Electrolysis Cells for High Yield Hydrogen Gas Production from Organic Matter[END_REF] or microbial electrosynthesis [START_REF] Rabaey | Microbial electrosynthesis-revisiting the electrical route for microbial production[END_REF]. Wastewater degradation in BES, like in anaerobic digestion, follows a multistep process: it includes biochemical disintegration of complex particulate matter to carbohydrates, proteins and lipids, followed by hydrolysis to sugars and amino acids and their fermentation to volatile fatty acids (VFA), carbon dioxide and hydrogen [START_REF] Velasquez-Orta | Evaluation of hydrolysis and fermentation rates in microbial fuel cells[END_REF].

Anode respiring bacteria (ARB) have shown ability to oxidize a wide range of substrates [START_REF] Pant | A review of the substrates used in microbial fuel cells (MFCs) for sustainable energy production[END_REF]. Acetogenic fermentation is usually the main metabolic pathway due to the higher ATP yield. Nevertheless, by increasing loading rate (more than 3 kg COD m -3 d -1 ) or reducing anode potential (less than -0.1mV vs NHE), this reaction can be inhibited due to hydrogen accumulation and mixed acid fermentation occurs [START_REF] Rabaey | Bioelectrochemical systems: a new approach towards environmental and industrial biotechnology[END_REF]. Hydrogen can be consumed by hydrogenotrophic methanogens and hydrogenotrophic ARB [START_REF] Kimura | Acetate oxidation by syntrophic association between Geobacter sulfurreducens and a hydrogen-utilizing exoelectrogen[END_REF]. Therefore, organic acid oxidation in electroactive biofilm (EAB) results in both a syntrophic relationship (fermenters/ARB) and a competition for fermentation byproducts (ARB/methanogens). This special microbial ecology relation may lead to the spatial heterogeneity of species and a reduction of coulombic efficiency due to losses in fermentation steps and methane production. Description of local segregation of fermentative and electroactive microbial groups is therefore essential to optimize the system design and enhance its efficiency in current production. Even if important understanding has been achieved in extracellular electron transfer (EET), up-scaling of BES technology is still facing major bottlenecks and industrial applications remain a future promise (Babauta et al., 2012).

It is believed that the optimization of substrate conversion efficiency and operating conditions are essential to maximize the BES performance [START_REF] Pandey | Recent advances in the use of different substrates in microbial fuel cells toward wastewater treatment and simultaneous energy recovery[END_REF]. In general, when supported by sufficient experimental data, numerical modelling could be a good tool to represent BES and simulate what operating conditions could lead to an optimum performance.

Therefore, mathematical modelling is a tool to integrate and understand complex interactions resulting from the coupling between transport phenomena, chemical and microbial reactions and the electrochemical processes occurring in BES. Several numerical models have been proposed in the last decade to describe and simulate a wide range of BES characteristics such as: the electroactive biofilm growth resulted from coupling microbial electron donor metabolism with anode as electron acceptor and electron transfer by conduction or mediator molecules [START_REF] Kato Marcus | Conduction-based modeling of the biofilm anode of a microbial fuel cell[END_REF][START_REF] Picioreanu | Mathematical model for microbial fuel cells with anodic biofilms and anaerobic digestion[END_REF][START_REF] Korth | A framework for modeling electroactive microbial biofilms performing direct electron transfer[END_REF]; EET mechanism at molecular scale [START_REF] Snider | Longrange electron transport in Geobacter sulfurreducens biofilms is redox gradientdriven[END_REF][START_REF] Renslow | Modeling biofilms with dual extracellular electron transfer mechanisms[END_REF], the role of diverse operating conditions (anode potential, pH influence, substrate limitation, hydrodynamic flow) [START_REF] Picioreanu | Mathematical model for microbial fuel cells with anodic biofilms and anaerobic digestion[END_REF][START_REF] Torres | Proton transport inside the biofilm limits electrical current generation by anode-respiring bacteria[END_REF][START_REF] Merkey | The performance of a microbial fuel cell depends strongly on anode geometry: a multidimensional modeling study[END_REF] or full BES behavior [START_REF] Pinto | Optimizing energy productivity of microbial electrochemical cells[END_REF]. One example of multispecies model EAB was proposed by Picioreanu et al. [START_REF] Picioreanu | Model based evaluation of the effect of pH and electrode geometry on microbial fuel cell performance[END_REF] and it considers syntrophy between various biomass metabolic types, integrating at the same time the calculation of pH with 2D liquid flow and transport of multiple solutes. Nevertheless, due to numerical cost of individualbased modeling approach (IbM) coupled to finite element solvers of mass, quantum and charge balances, this approach requires large computing capacity and, thus, its application is limited for optimization procedure.

Glucose is a substrate of particular interest because it is the smallest monosaccharide to be fermented to VFA. The syntrophic glucose oxidation in microbial fuel cell has been experimentally studied since the 1980s [START_REF] Bennetto | Anodic reactions in microbial fuel cells[END_REF] and furthermore the metabolic degradation pathway is a well-described process [START_REF] Freguia | Syntrophic processes drive the conversion of glucose in microbial fuel cell anodes[END_REF]. A thermodynamic approach can be used to estimate conversion yields and maximum biomass specific rates, because it offers more consistency and generality when working with of different (unknown) microbial metabolic groups [START_REF] Heijnen | Bioenergetics of Microbial Growth[END_REF]. In order to estimate segregation of microbial populations in EAB, modelling trends can be compared with experimental results obtained from a cascadefed system [START_REF] Hodgson | Segregation of the Anodic Microbial Communities in a Microbial Fuel Cell Cascade[END_REF].

The aim of this work was to develop a 2D multiphysics mechanistic numerical modeling respecting the thermodynamic laws of microbial metabolism in order to simulate at cellular scale, inter-species segregation in mixed anodic biofilms. The proposed approach relied on the assembly of the afore-mentioned numerical approaches in order to put forward a global description of BES at cellular scale. For this purpose, one particularity of the numerical model developed was to make a spatiotemporal variation of the biofilm growth zone by using a moving mesh according to the arbitrary Lagrangian-Eulerian method (ALE). Segregation have been simulated both in the biofilm depth and along the biofilm length, when the biofilm is fed with a solution of glucose flowing over the biofilm surface, under various operating conditions.

Furthermore, the influence of several operational parameters (i.e., inlet concentration, poised anode potential, specific surface area and retention time) on the biofilm microbial composition has been investigated. Finally, thanks to modelling, computed coulombic efficiency and organic removal rate allowed to define an optimal Pareto front helping to design and to select best operating conditions for bioanodes. This multi-objective optimization method is still few used in BES development.

Model description

Overview

The usual assumption is to consider that the most important gradients (concentrations, pH, electrical potential) are in the direction from bulk liquid to the support material (i.e., anode), thereby constructing one-dimensional (1D) biofilm representations. Sometimes, this assumption has to be relaxed, for example in case of an anode with irregular geometry (porous electrode, [START_REF] Merkey | The performance of a microbial fuel cell depends strongly on anode geometry: a multidimensional modeling study[END_REF][START_REF] Picioreanu | Model based evaluation of the effect of pH and electrode geometry on microbial fuel cell performance[END_REF]) or when microbial interactions at colony level are studied (interspecies electron and substrate transfer, [START_REF] Picioreanu | Mathematical model for microbial fuel cells with anodic biofilms and anaerobic digestion[END_REF][START_REF] Picioreanu | Model based evaluation of the effect of pH and electrode geometry on microbial fuel cell performance[END_REF]. Several model parameters have been identified as essential for the representation of experimental data: biofilm thickness, biomass concentration, electron transfer properties (biofilm conductivity, cofactor redox concentration, mediator diffusivity) [START_REF] Kato Marcus | Conduction-based modeling of the biofilm anode of a microbial fuel cell[END_REF][START_REF] Korth | A framework for modeling electroactive microbial biofilms performing direct electron transfer[END_REF], effective diffusion of diverse solutes in the biofilm (which influence substrate depletion and local acidification) [START_REF] Renslow | Modeling biofilms with dual extracellular electron transfer mechanisms[END_REF][START_REF] Torres | Proton transport inside the biofilm limits electrical current generation by anode-respiring bacteria[END_REF], metabolic activity (maximum biomass growth rate, maintenance rate, EPS production) [START_REF] Merkey | The performance of a microbial fuel cell depends strongly on anode geometry: a multidimensional modeling study[END_REF][START_REF] Jayasinghe | Metabolic modeling of spatial heterogeneity of biofilms in microbial fuel cells reveals substrate limitations in electrical current generation[END_REF]. However, the local segregation of diverse microbial metabolic groups in EAB due to interspecies competition (electrogens, acidogens, acetogens, methanogens, etc.) has been less approached in numerical studies. The present numerical approach strives for conversion step model of an organic compound into electrons and the microbial interspecies competition in a biofilm growing along an anode over which a solution of glucose flows. For clarity reasons, operating conditions are considered as isothermal (298K) even if the biochemical growth model used is able to represent temperature variations using correction factor [START_REF] Heijnen | Bioenergetics of Microbial Growth[END_REF][START_REF] Gildemyn | A Gibbs Free Energy-Based Assessment of Microbial Electrocatalysis[END_REF]. For this scope, a two-dimensional (2D) rectangular domain of fixed dimensions was constructed. The computational domain was further divided into two regions of variable geometry (Figure 1): a biofilm area (BF) and a bulk liquid area (BK), spanning over the anode length (direction x) with a variable thickness (direction y). While convection and diffusion of solutes takes place in the bulk liquid (due to flow over the biofilm), diffusion, reactions and electron transfer occur in the biofilm. The essential solutes considered in this model are: glucose, acetate, hydrogen and methane. The biofilm consists of five biomass types: glucose consumers (fermenters), hydrogen consumers to methane (methanogens), electroactive acetate and hydrogen utilizers (acetotrophic and hydrogenotrophic electroactive bacteria) and inert biomass resulted from inactivation. The general pathway of glucose conversion follows this framework:

1) Mass transport of solutes in the bulk liquid flowing over the biofilm (diffusion and convection)

2) Mass transfer of solutes to/from the biofilm through the biofilm/liquid interface (flux continuity)

3) Conversion processes in the biofilm volume. These include: acidogenesis (glucose consumption with acetate and hydrogen production resulting in biomass growth), methanogenesis (growth on acetate consumption to methane and CO 2 )

and electroactive acetate and hydrogen utilization by ARB.

4) Electron transfer to the anode via an equivalent conductive biofilm matrix.

When current production is the process goal, the methane production by hydrogenotrophic microorganisms is an unwanted side-reaction that consumes hydrogen and leads to a reduction of coulombic efficiency. Following the EAB modelling approach previously proposed [START_REF] Richter | Cyclic voltammetry of biofilms of wild type and mutant Geobacter sulfurreducens on fuel cell anodes indicates possible roles of OmcB, OmcZ, type IV pili, and protons in extracellular electron transfer[END_REF], the process by which the glucose feed is converted into electrons can be split into several steps (Figure 1A):

1) glucose feed and transport in liquid; 2) glucose transfer to biofilm; 3) glucose transport in biofilm; 4) glucose fermentation; 5) acetate and hydrogen transport in biofilm; 6) acetate and hydrogen conversion by methanogens and electrogens; 7) electron transfer to anode. This study does not aim to fit biofilm intrinsic parameters, but rather to propose a conceptual framework for optimization of a biofilm system with segregated biomass activity. This is why the metabolism of different microbial groups is described using the Gibbs energy approach [START_REF] Heijnen | Bioenergetics of Microbial Growth[END_REF]) and a compilation of other parameters from the literature. Even if a bioenergetic theory applied to BES has been proposed elsewhere [START_REF] Gildemyn | A Gibbs Free Energy-Based Assessment of Microbial Electrocatalysis[END_REF], it is clear that some of the parameters implemented (i.e., maintenance, biofilm conductivity or fermenter hydrogen inhibition) are highly sensitive and can lead to strong variations in results. The parameters values used in this study are selected with regards to the literature review given in Figure S4 (supplementary material). In particular, we took into consideration the highly reproducible results obtained with rather similar conditions by researchers from the University of Gent [START_REF] Zhang | Rapid and quantitative assessment of redox conduction across electroactive biofilms by using double step chronoamperometry[END_REF]. Operational parameters have been chosen to take only into account the influence of reactor design (anode length, bulk liquid height, specific surface) and the conditions monitored by the operator (hydraulic retention time, anode potential, glucose concentration). Indeed, this work aims to propose a framework for BES operational parameters optimization and not a biofilm intrinsic parameters study.

Nevertheless, under constant microbiological description (Table 1A), the model can be used to test various experimental operating conditions (Table 1B) and can suggest general trends of local segregation and biofilm efficiency.

Biofilm

BF is a 2D component (Figure 1B-C). It describes biofilm growth, and thus, interspecies competition, local potential and substrates conversion and consumption.

It is geometrically defined by the anode length (L x ) (x-axis) and the biofilm thickness (L f ). In order to represent biofilm growth, L f is modeled using a moving mesh based on an Arbitrary Lagrangian-Eulerian method. This allows a deformation of the meshing under constraints while mesh number stays constant (Nx,f=100; Ny,f=50).

Substrates conversion model

Syntrophic glucose conversion [START_REF] Freguia | Syntrophic processes drive the conversion of glucose in microbial fuel cell anodes[END_REF] is used to describe main chain reaction in biofilm (Figure 1A). Table 2 describes the different metabolism steps used in the model to represent glucose oxidation. For each step, a detail of catabolic and anabolic reactions is calculated by using bioenergetics growth theory. In order to simplify model, all biomasses are considered identical with the same density and same composition formula (CH 1.8 O 0.5 N 0.2 ).

Table 2 Catabolic reaction produces energy to run biomass production (anabolic reaction). However, a part of the energy is dissipated which depends mainly on substrate characteristics (carbon chain length and oxidation number) [START_REF] Heijnen | A black box mathematical model to calculate auto-and heterotrophic biomass yields based on Gibbs energy dissipation[END_REF]. Catabolic and anabolic energies are function of temperature and substrates and byproducts concentrations. Considering catabolic reaction as an energy pump with dissociation losses, solutes conversion yield can be calculated for each term of each step (Y i,Xi ). Details of this calculation are provided in Supplementary Material.

Growth kinetic model

Kinetics used to describe biomass growth rate are all based on Monod law. In fermenters metabolism, H 2 partial pressure can inhibited dehydrogenase enzyme responsible of NADH [START_REF] Fukuzaki | Inhibition of the Fermentation of Propionate to Methane by Hydrogen, Acetate, and Propionate[END_REF]). An inhibition factor in glucose fermentation step is introduced in an extension of Monod equation taking into account hydrogen concentration (eq. 4) [START_REF] Han | Extended monod kinetics for substrate, product, and cell inhibition[END_REF]:

μ /0 = μ 12+,/ 1 - 3 3 4 5 $ 6 + 7 ,/ 4
C H2 inhibit is the concentration of hydrogen above which reaction is inhibited (in mmol.L - 1 ), assumed equal to the solubility in water at 298K (0.8mmol.L -1 ) [START_REF] Lide | CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data[END_REF].

Methanogens growth is taken from a simple Monod form depending on hydrogen concentration (eq. 5):

μ /1 = μ 12+,/1 3 3 + 7 3 ,/1 5
For ARB, growth rate are expressed using a Nernst-Monod law [START_REF] Kato Marcus | Conduction-based modeling of the biofilm anode of a microbial fuel cell[END_REF], which depends on electron donor concentration (acetate (X a ) or hydrogen (X h ))

and local acceptor potential (eq. 6):

μ / = μ 12+,/ + 7 ,/ 1 1 + :; - < = -> ?2,/ @ 6
E ka corresponds to the potential where the current is half of the maximum current in a typical sigmoid shape cyclic voltammetry. It represents the EET ability and can be related to the formal potential of the redox species X i . This is true for a pure culture grow in specific condition. Authors are aware that EET is dependent on conditions growth such as anode potential or substrates [START_REF] Levar | An Inner Membrane Cytochrome Required Only for Reduction of High Redox Potential Extracellular Electron Acceptors[END_REF] and that in mixed culture several redox systems can operate at the same time [START_REF] Rimboud | Multi-system Nernst-Michaelis-Menten model applied to bioanodes formed from sewage sludge[END_REF]. However, in order to limit model complexity, two constant EET will be considered (represented by two E ka values taken from the literature [START_REF] Zhang | Rapid and quantitative assessment of redox conduction across electroactive biofilms by using double step chronoamperometry[END_REF] in pure culture and associated to one of each electroactive biomass).

From equations 4-5, we can deduce solutes conversion rate in each step (eq. 7) and electron conversion (eq. 8):

& ,/ = B ,/ C / D μ / = B ,/ E / F D μ / 7 & ,/ = B ,/ C / D μ / = B ,/ E / F D μ / 8
With r i, Xi the electron donor i consumption rate for biomass X i (mol m -3 s -1 ), c X the biomass molar density (mol m -3 ), Y i,Xi the stoichiometric conversion yield solutes i in biomass X i (mol i mol X -1 ), ρ X the biomass density (g m -3 ) and M (g mol -1 ) the biomass molar weight. Introducing the hypothesis of a constant limitation in intracellular electron transport chain in microorganisms (cf. Supplementary Material). Heijnen et al.

proposed a method to calculate the maximum consumption rate and maximum biomass specific growth for each reaction [START_REF] Heijnen | Bioenergetics of Microbial Growth[END_REF]. These parameters depend on the catabolic reaction energy and the maintenance need and, thus, on temperature and concentration. Even if limitations of this assumption can be discussed, it allows a reduction of model degrees-of-freedom number and avoid a supplementary correction factor. No growth can be observed if catabolic energy production is lower than the maintenance energy. Therefore, for each reaction, a threshold substrate concentration to sustain microbial growth can be determined for fermenters and methanogens (eq. 9), and, similarly, for the electroactive bacteria, including local potential dependence (eq. 10).

With m i , the specific rate of electron donor is catabolized to generate the necessary Gibbs energy to maintain the biofilm (in mol i s -1 ). Below these values, biomass will not be able to grow and auto-consumption maintenance will be forced (cf. biomass balance section).

Resource competition model [START_REF] Grover | Resource Competition[END_REF] proposes a useful framework to predict composition of a multispecies population competing for the same substrate. In long term study, the species having the lowest threshold concentration will be dominant. In this study, this approach has been applied between electroactive hydrogenotrophic bacteria and hydrogenotrophic methanogens. As hydrogen level concentration remains low in the biofilm and no resource is sequestered in bacteria, the less limited species will eventually dominate the other. Therefore, a simplification hypothesis is introduced to force only the growth of the most favorable biomass (lower C min ) while the other is totally inhibited.

Solutes mass balance equations

In order to reduce model complexity, all the components are considered as part of the liquid phase. Solutes concentration distribution is represented by the Fick law (eq. 11). Indeed, no convection within the biofilm is supposed and electro-migration is neglected as bulk is considered as well-buffered. Diffusion takes place in both x and y axis.

]) * ]6 = ∇ . " 2^6 " ∇ + & = ∇ . " 2^6 " ∇ + ∑ & ,/` (11)
At the biofilm/bulk interface, concentrations continuities are ensured allowing exchange from the two domains. However, as solutes mobility is reduced by biofilm structure, a correction factor (D fact ) of 0.8 was applied to the diffusion coefficient in water inside the biofilm [START_REF] Kissel | Numerical Simulation of Mixed-Culture Biofilm[END_REF]. At the electrode/biofilm interface, an insulation boundary condition was applied as no solutes can flow through the solid electrode (

Charge balance equations

As electron balance time characteristic is about 1000 times smaller as microbial growth [START_REF] Kissel | Numerical Simulation of Mixed-Culture Biofilm[END_REF], a pseudo-steady-state was assumed and a Poisson law was used to represent it (eq. 12).

0 = ∇ . T σ efg ∇ V -E jk Y + r Q,mk + r Q,mn (12)
With σ bio the biofilm conductivity (S m -2 ), r e,Xa the electron source term in Xa (A m -3 )

and r e,Xh the electron source term in Xh (in A m -3 ). The source term r e represents the local current production by ARB in the biofilm. It depends on conversion rates (Y e, Xi ), and thus to catabolic energy available (cf SI). Biofilm conductivity was considered as constant in the biofilm, representing a fully conductive matrix with a transfer only orientated in y-axis [START_REF] Kato Marcus | Conduction-based modeling of the biofilm anode of a microbial fuel cell[END_REF]. Anode was considered as an equipotential surface. Thus, at anode/biofilm interface, a Dirichlet boundary condition is considered to represent the anode poised potential (V(y=0) = V an ). Anode was supposed to be perfectly conductor; insulation conditions were applied of the biofilm matrix ( = 0 . Then, we can calculate, the local current density collected at the anode/biofilm interface using the ohm law (eq. 13), the total current production in the bioanode (eq 14) and the mean current density in the bioanode (eq. 15).

Biomass balance equations

Biofilm is composed of 5 different volume fractions:

-4 actives (X g , X a , X m , X h ), where the four conversion steps in glucose syntrophy take place, -1 inactive (X in ), representing dead cell, EPS fraction and insoluble byproduct.

Biomass composition is based on local biomass volume fraction using "fuzzy layer" approach [START_REF] Wanner | A multispecies biofilm model[END_REF]. At any times, anywhere in the biofilm, the sum of all biomass volume fractions is constant and validates (eq. 16):

~D = 1 16

Biomass growth is limited to one direction (y-axis) perpendicular to electrode surface.

A convection diffusion equation was used to represent biomass balance in biofilm (eqs. 17-18). As defined earlier, biomass activity can be separated in two domains modelled by a Boolean condition. If maintenance needs are met (C i > C min Xi ), biomass can grow and invest energy in EPS matrix production. Otherwise, active biomass will be progressively converted into inactive biomass, reducing the biofilm ability to oxidize substrate and/or produce electron.

D ! -∇ . " / ∇ D + •D € = & / = • 1 -' ƒ"… C / D μ / †' > ‡fˆ/ -‰ 42 C / D †' < ‡fˆ/ 17 D 4 ! -∇ . " / ∇ D 4 + •D € = & / 4 = ~•' ƒ"… C / D μ / †' > ‡fˆ/ ‰ 42 C / D †' < ‡fˆ/ Š / -‰ a‹ D 4 18
With D X the biomass diffusion coefficient (voluntary low) just to help system solution convergence (mol X m -2 s -1 ), v the local advection in biofilm (m s -1 ), f EPS the biomass growth fraction invested in EPS production, M the biomass molecular weight. In this study, inactive biomass lysis was not considering (b lys =0). Growth of the inner layers of the biofilm creates an advection movement. Therefore, we can calculate the local advection velocity ν as the sum of the contribution of each volume fraction in the inner layers using eq. 19:

OE = 1 - x x 12+ • & /0 + -Ž b & /2 + & /5 + & /1 + & / 4 {€ 19 Using the coefficient •1 - r Ž -ŽLMN
', a maximum biofilm thickness is imposed regardless of the physical limiting factor (shear stress detachment, Red/Ox concentration gradient saturation). Advection velocity at the outside bulk/biofilm interface (ν(L f )) was applied as a moving mesh boundary condition while biomass boundary condition was considered as an outflow. Therefore, the detachment only concerns the outermost layer of the biofilm.

Kinetics reaction for solutes concentration, biomass fraction and potential inside the biofilm are summarized in table 3.

Table 3 2.4. Solving method

The model has been set up with COMSOL Multiphysics® 5.2. Biofilm has been divided in 100 meshes along x-axis and 40 along y-axis. X-meshes were constant along anode whereas Y-meshes had a symmetric distribution with a 5-factor size between center and boundaries in order to get best resolution at the interfaces. A progressive boundary refinement was considered to represent the bulk/biofilm interface.

Simulations were achieved using a PARDISO direct solver computing a system of linear equation with 33720 degrees of freedom representing the total set of equations. The study was divided in 2 steps. First step was a stationary study that solves solutes concentration (BK, BF) and local potential (BF) and uses a highly non-linear Newton Raphson solver method. This allows to reach better initial conditions. Second step was a time dependent study incorporating biomass growth. The model was solved with a fully coupled method using an automated damping factor adjustment Newton Raphson method. Time dependent solver use a backward differentiation method. Time steps were taken freely by the solver using a non-linear controller in order to avoid oscillations and help numerical convergence (detailed solving method chronology are available in figure 2). Physical time of fifteen days simulation was usually enough to reach a biological steady state in most operating conditions which is relevant with literature [START_REF] Ren | Characterization of Microbial Fuel Cells at Microbially and Electrochemically Meaningful Time scales[END_REF]. Simulation solving times varied from 0.3 to 3h depending on the tested conditions. Fig. 2.

Results and discussion

Main trend of simulated biomass behavior (model validation)

The model is able to simulate biofilm growth and composition along anode length. Figure 3 shows an overview of the main model results of a simulation obtained with a 5 cm anode length and 3mmol.L -1 glucose inlet concentration. A one-hour Hydraulic Retention Time was imposed (u flow = 0.13 mm.s -1 ) to highlight species segregation.

Fig. 3.

Bulk concentrations are indicators of solute accumulation in the biofilm (Figure 3A). If bulk solute concentration increases, it means that byproduct excess is transfer to the bulk. In contrast, when bulk concentration decreases, solutes are more consumed than produced in the biofilm. Glucose concentration decreased along anode length. Acetate only accumulated at the entrance and then it was all consumed inside the biofilm (c b,ac =0). Hydrogen accumulated in the first centimeter but was also consumed along the anode length. Glucose fermenter were mainly concentrated in the outer layer of the biofilm where glucose concentration was higher. In contrast, ARB (acetate-fed and hydrogen-fed) were concentrated in the inner layer due to local potential limitation. Indeed, even if anode poised potential was favorable to run metabolism (E ka < V an ), potential ohmic dropped due to biofilm conductivity limits ARB growth to 30-40µm. Methanogens were almost not represented in this case because conditions were favorable for hydrogen ARB (high potential) which consumed almost all the hydrogen available. In addition, as outermost layer suffers from detachment (biomass outflow) no ecological niche exists for methanogens. In this study, hydrogenotrophic ARB Monod constant have been chosen arbitrarily low which leads to a faster growth of this biomass fraction and a predominance in the inner layer.

To deal with the whole issue of geometric parameter and operating conditions impact on biomass behavior, main parameters such as anode potential, inlet flow rate, glucose concentration, channel size and hydrogen threshold have been investigated.

Anode potential variation study is reported in Supplementary Material. Therefore, the model was able to represent switches of different electroactive bacteria activities related to anode potential from a reference condition. It is important to note that these switches were directly linked to intrinsic kinetics parameters (mid-term potential, Monod constant) which are quite stable and reproducible for electroactive bacteria in pure culture like Geobacter under well-controlled conditions but need to be assessed more precisely in mixed culture fermenting bioanode. In addition, higher flow and higher channel length lead to a higher glucose conversion (no inhibition) and a higher fermenter volume fraction in the biofilm competing with electroactive volume fraction and thus reducing current production. Contrary, higher initial concentration lead to a faster inhibition due to hydrogen concentration (cf. Supplementary Material).

The over hydrogen threshold glucose conversion into electrons was limited by fermenting rate. Indeed, despite of favorable hydrogenotrophic metabolism rate (low Monod constant, high anode potential), hydrogen has accumulated and inhibited glucose fermenter metabolism. As a first approximation, hydrogen inhibiting concentration has been evaluated as hydrogen solubility in water [START_REF] Lide | CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data[END_REF].

Influence of hydrogen threshold on glucose conversion and electroactivity is reported as additional materials (cf. Supplementary Material). Increasing hydrogen threshold has accelerated glucose conversion and increased fermenter part in biofilm volume.

Direct comparison between our mathematical model and experiment reported in the literature was also achieved. Freguia et al. studied glucose fermentation in a microbial fuel cell with graphite bed electrode [START_REF] Freguia | Syntrophic processes drive the conversion of glucose in microbial fuel cell anodes[END_REF]. They tested the influence of poised anode potential on cell mass balance and current production, which gave information on syntrophic behavior in anodic biofilm. In our model, the graphite bed anode has been considered as a 15 cm length flat plate with a 10mm channel height (Ly). Under model hypothesis, only 35% of glucose was converted in acetate and hydrogen after a ten days' simulation which was likely due to three reasons i) the limited specific surface area to comply convective flow conditions ii) the absence of suspended cell activity iii) the hydrogen concentration inhibition. A low acetate and hydrogen concentrations were present in the outlet which was relevant with experimental data. However, simulated methane outlet was a two order higher than observed one. This could be due to the presence of another hydrogen consumption pathway or an inhibition of hydrogenotrophic methanogens by back-diffusion of oxygen in real conditions. Current production simulated was two order smaller than measured which is directly linked to the low glucose conversion rate. In conclusion, lack of information and control on real operating conditions -compared to ideal model approach-lead to variation on absolute values of performance indicators. However, metabolism simulated is relevant with experimental data concerning conversion pathways general trends.

Cascade-fed MFC system gives details of the segregation along anode length [START_REF] Hodgson | Segregation of the Anodic Microbial Communities in a Microbial Fuel Cell Cascade[END_REF]. Hodgson et al presented a system of four MFCs in series, with a constant external resistance of 20Ω. Feedstock was a dry dried distillers' grain medium compound of glucose, pentose, glycerol and lactate with a total carbohydrate of 16.8 g L -1 (glucose equivalent). Solutes concentrations, current density and microbial community analysis have been reported. In the model, only one fermenting pathway was considered, therefore fermentable carbohydrates byproduct (except acetate)

were considered as a glucose equivalent concentration (on a C_mole basis). Four simulations have been run successively with previous bulk outlets concentrations as next reactor inlets. In these simulations, no inhibiting fermentation factor has been considered. Indeed, preliminary tests (not reported) have shown a 5% glucose removal with a 0.8mM inhibiting concentration which is much lower than observed one.

Current production in model has increased, like in experimental study, along the cascade MFC. However, it was a two order higher with similar specific surface area considered. This is likely due to two model hypothesis i) As exposed before, no recombination and constant and perfect transfer was supposed through the whole biofilm thickness ii) the kinetic model of ARB has been defined from Geobacter studies.

Indeed, Geobacter are able to develop thicker active biofilm and perform long range extracellular electron transfer.

From these direct comparisons, relevant evolution of biomass fraction was observed compare to experimental data: fermenting bacteria volume fraction has decreased while electroactive bacteria volume fraction increased. Therefore, syntrophy between various biomass metabolic types was enough well modelled in order to run "numerical" experiences to attain optimal coulombic efficiency.

Toward multi-objective optimization

The main interest of modelling approach appears when process behaviors can be predicted and more efficient design proposed. The final aim of this work, using previously determined metabolism properties, is to propose a framework to set up multicriteria studies for design optimization. For an efficient bioanode reactor, main objectives to increase performance would be to maximize flow and conversion from glucose to electron, i.e. maximize coulombic efficiency and organic removal rate. Using the model, influence of several cross-parameters can be tested and gives information on reactor performance. Four parameters were varied to check their individual sensitivities on current densities and solutes conversion (cf. Supplementary Material).

Coulombic efficiencies and organic removal rates were represented as multi-objective optimization (Figure 4). Low channel height should be favored to increase both coulombic efficiency and organic removal rate. Similar trends have been observed with increasing glucose inlet concentration and anode potential. However, threshold influence has to be determined. On the one hand, higher anode potential mainly has increased the coulombic efficiency. On the other hand, higher glucose inlet concentration mainly has increased the organic removal rate. Furthermore, high organic removal rate was reached for low HRT, however, the coulombic efficiency could be increased by increasing HRT. The anode length and bulk velocity (linked to HRT) have, thus, to be optimized with regards to the main objective aimed.

According modelling approach and experimental results, it is possible to define the following hierarchy: i) Anode potential (V an )

ii) Operating conditions inlet inflow (u flow ) and concentration of glucose (C Gl 0 )

iii) Geometry parameters channel height (Ly) and anode length (Lx)

Fig. 4.

Anode potential is the main factor influencing coulombic efficiency acting on electroactive biofilm growth. This parameter has been considered at a constant value of 0.2V in a multi-criteria analysis in order to do a screening of operating parameters where a Pareto front brings in light best design conditions (Figure 5). Under simulated conditions, highest coulombic efficiency was closed to 65-70% which is relevant with glucose-fed bioanode experiment [START_REF] Ichihashi | High-Performance Bioanode Development for Fermentable Substrates via Controlled Electroactive Biofilm Growth[END_REF]. However, prediction accuracy of such approach is directly linked to the number of simulations. Extensive screening of operating conditions should be run to attest parameter interdependence and would be presented in future work. This model framework could be extended in order to propose new design to optimize flow along reactor in various type of BES. As an example, it could be applied to improve performance in cascade fed reactors for electrons conversion (MFC, MEC) or byproducts synthesis (electro-fermentation).

Conclusion

This work proposes a 2D multispecies glucose-fed bioanode model working under continuous flow. Using glucose syntrophy hypothesis, the model is able to represent biomass segregation, solutes concentration distribution from the inlet to the outlet and current production. Influence of hydrogen inhibiting concentration on glucose fermentation clearly appeared as a limiting factor for overall conversion (glucose to electron).

We propose a framework for BES design optimization using a multicriteria approach (coulombic efficiency VS organic removal rate). Under a set of hypotheses for metabolism description and kinetics, multi objectives Pareto front was exhibited in order to determine optimal design and operating conditions. Conceptual model for the segregation of microbial metabolic types within the biofilm depth and along the flow direction (inert biomass is not represented here). C) Mathematical model set up with bulk liquid and biofilm sub-domains. There is a laminar flow in the liquid bulk (u flow ) supporting the transport of solutes. The variable biofilm thickness (L f ) results from a moving mesh driven by microbial growth inside the biofilm with different rates based on substrate and electron transfer availability along the anode length. Fig. 4: Coulombic efficiency and organic removal rate for several operating conditions (Blue: Hydraulic Retention Time (HRT = 0.1, 0.3, 0.5, 0.8, 1h) ; Orange : glucose inlet concentration (C i Gl 0 = 1, 2, 3, 4, 5mM) ; Green : channel height (L y = 1.5, 1.9, 3 ,5, 7, 10mm); Yellow: Anode potential (V an =-0.2, 0, 0.1 0.2V); Red: reference: HRT: 1h, C i Gl 0 : 3mM; L y : 1.9mm, V an : 0.2V)).

Fig. 5: Multi-criteria analysis for multi-objectives optimization for several operating conditions (HRT = 0.1, 0.5 and 1h ; c i Gl 0 = 1, 2, 3, 4 and 5mM ; L y = 1.5, 1.9 and 3 mm

Fig

  Fig. 1.

  BK is a 2D component which represents the feedstock influent flux along the biofilm/liquid interface (Figure1B-C). The channel is a framework defined by the anode length (Lx) and the bulk height (Ly). Inlet is supposed to have an average velocity constant (u flow ) with sufficient entrance length in order to ensure a stable fluid velocity distribution in x=0. A laminar flow is ensured in flow speed and bulk height ranges calculating dimensionless Reynolds number < 100 . Convective flow is maintained superior to diffusive flow calculating dimensionless Péclet number > 1 , thus, there is no solutes accumulation in bulk. These conditions permit to reduce model complexity and thus, meshing detailed. In the model, 100 meshes along anode length (x-axis) and 5 meshes along bulk height (y-axis) are considered. However, pure diffusive flow could be addressed with a consequent extra numerical calculation demand (especially for long anode length). Fluid is considered as incompressible and behaves like a Newtonian fluid with a constant viscosity. A simplified Navier-Stokes equation (eq.1) and global mass balance (eq.2) can be used to model flow rate[START_REF] Duddu | A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment[END_REF].Solutes concentrations are calculated via a convection/diffusion mass balance equation (eq. 3):! = ∇ . " , ∇ ,$ + . ∇ ,$ -& ,$3Inlet (x=0)is only composed of glucose (except for cascade fed set up) and is maintained constant during each simulation. The outer boundary layer (x=L x ) is an outflow ( " () * (+ +,-+ = 0 ). In order to reduce model complexity, influence of suspended cell bacteria has been neglected (r i,b = 0). Thus, solutes conversion only takes place inside the biofilm. Continuities of flux and concentration are considered at bulk/biofilm interface.

  were also supposed on the vertical backside of the biofilm.

Fig. 5 .

 5 Fig. 5.

Fig. 2 .

 2 Fig. 2. Solving method scheme.Fig. 3. A) Solutes concentrations at bulk/biofilm interface along anode length (x-axis) . (red line: glucose; green line: methane; blue line: acetate, cyan line : hydrogen). B) Biofilm cross sections (in x=0.5mm, x=5mm, x=25mm, x=45mm) showing solutes concentrations (red line: glucose; green line: methane; blue line: acetate, cyan line: hydrogen) and biomass volume fraction (red line Xg; blue line: Xa; green line: Xm; black line: Xin).

Table 1 2

 1

.2. Bulk

Table 1 :

 1 List of model parameters. A) Intrinsic biofilm properties (kept constant during the simulations). B) Monitoring external parameters (can be varied as model input).

Table 2 :

 2 Catabolic et anabolic reactions in biofilm from[START_REF] Heijnen | Bioenergetics of Microbial Growth[END_REF]

Table 3 :

 3 Stoichiometric matrix of kinetics reaction in the biofilm Fig.1. Schematic model description. A) Species syntrophic relation and transport chain with reactions. Xg: glucose consumers; Xm: methanogens; Xh and Xa: electrogens. B)

  ). Insert: Conceptual multi-objective Pareto front representation.

	A	Parameter description Step Reaction	Symbol Value Stoichiometric reaction Units	Source
	Biofilm properties Initial biofilm thickness Maximal biofilm thickness Glucose fermentation by fermenter (Xg) Molecular weight of a C-mol biomass Total dry biomass weight density in the biofilm Catabolic Anabolic Total biomass molar density in the biofilm Biofilm conductivity Catabolic Acetate consumption by ARB (Xa) Anabolic	Lf init Lf max Μ ρX CX σbio 0.525	10 60 24 50 2032 + 4 0.175 → 0.4 0.001 + 4 → 2 + 0.2 µm µm g mol -1 + . . g L -1 mmol L -1 mS cm -1 → 2 + 0.2 + 0.275 Based on an elemental composition of Assumed Compared in S4 + 2 + 4 + 4 + 0.25 . + 0.05 biomass CH1.8O0.5N0.2 (Korth, et al. 2015) Calculated Compared in S4 + 9 + 8 → . . . + 0.05 + 0.4
	Microbial kinetics and thermodynamics Inactivation rate EPS production coefficient For fermenter Xg Catabolic Hydrogen consumption by ARB (Xh) Anabolic	bina feps	0.1 0.1 2.1 + 0.2	day -1 -→ 2 + 0.8 →	+ 2 + . .	Assumed Assumed . + 2.5
	Anabolic Gibbs energy a Catabolic Gibbs energy a Dissipation Gibbs energy a Monod constant Hydrogen consumption by methanogens (Xm)	Catabolic Anabolic	∆G 01 ∆G 01 ∆G 01 diss,Xg ana,Xg cat,Xg KGl Xg	-245 -206 -236 2.6 + 0.25 2.1 + 0.2 kJ molXg -1 kJ molGl -1 kJ molXg -1 mol m -3 + 0.25 + 0.8 → 0.25 (Heijnen and Kleerebezem 2010) Calculated Calculated (Batstone, et al., 2002) + 0.75 + → . . . + 2.5
	Inhibiting hydrogen concentration		CH2 inhibit	0.8	mmol L -1	Discussed in S3
	For ARB Xa
	Anabolic Gibbs energy a		∆G 01			ana,Xa	29	kJ molXa -1	Calculated
	Catabolic Gibbs energy a		∆G 01	cat,Xa	-137	kJ molAc -1	Calculated
	Dissipation Gibbs energy a		∆G 01	diss,Xa	-421	kJ molXa -1	Calculated
	Monod constant		KAc Xa	0.1	mol m -3	Compared in S4
	Mid-term potential		Eka Xa	-0.1	V (Vs SHE)	Assumed
	For ARB Xh
	Anabolic Gibbs energy a		∆G 01			ana,Xh	25	kJ molXh -1	Calculated
	Catabolic Gibbs energy a		∆G 01	cat,Xh	-102	kJ mol H2	-1	Calculated
	Dissipation Gibbs energy a		∆G 01	diss,Xh	-986	kJ molXh -1	Calculated
	Monod constant		KH2 Xh	0.008	mol m -3	Assumed
	Mid-term potential		Eka Xh	0.116	V (Vs SHE)	(Kimura and Okabe, 2013)
	For methanogen Xm
	Anabolic Gibbs energy a		∆G 01		ana,Xm	25	kJ molXm -1	Calculated
	Catabolic Gibbs energy a		∆G 01	cat,Xm	-34	kJ molH2 -1	Calculated
	Dissipation Gibbs energy a		∆G 01 diss,Xm	-986	kJ molXm -1	Calculated
	Monod constant		KH2 Xm	0.0008	mol m -3	(Batstone, et al., 2002)
	Mass transfer coefficient
	Glucose diffusion coefficient in water		DGl	0.5⋅10 -9	m² s -1	(Batstone, et al., 2002)
	Bulk acetate diffusion coefficient		Dac	1.2⋅10 -9	m² s -1	(Batstone, et al., 2002)
	Bulk hydrogen diffusion coefficient		DH2	5.0⋅10 -9	m² s -1	(Batstone, et al., 2002)
	Bulk hydrogen diffusion coefficient		DCH4	1.5⋅10 -9	m² s -1	(Batstone, et al., 2002)
	Dynamic viscosity		µdyn	0,896⋅10 -4	kg m -1 s -1	(Lide, 2001)
	Faraday constant		F	96485	C mol -1	(Lide, 2001)
	Gas constant		R	8.314	J mol -1 K -1	(Lide, 2001)
	B	Parameter description					Symbol	Range	Units
	Cell design
	Anode Length						LX	1-10	cm
	Bulk liquid height						LY	0.05-2	cm
	Liquid flow section area						Ab	0.05-50	cm²
	Bulk liquid volume						Vb	0.05-750	cm 3
	Specific surface area (1/Ly)						Sa	2500	cm -1
	Operating conditions
	Temperature						T	298	K
	Hydraulic retention time						HRT	0.5 -10	h
	Poised anode potential						Van	-0.2 -0.2	V
	Inlet glucose concentration						CGl 0	1 -10	mmol L -1
	Inlet acetate concentration						CAc 0	1.10 -5	mmol L -1
	Inlet hydrogen concentration						CH2 0	1.10 -5	mmol L -1
	Inlet methane concentration						CCH4 0	1.10 -5	mmol L -1

a Gibbs energy are considered in the biochemical standard state

(25 °C, pH 7, other concentrations 1M)

Acknowledgements

This research was partially supported by a PhD fellowship from the French ministry of Research.

Authors gratefully thank Dr Cristian Picioreanu from TU-Delft for assistance in numerical modelling with Comsol, and for comments that greatly improved the manuscript.

Appendix A. Supplementary

Supplementary material of this work can be found in online version of the paper References Babauta, J., Renslow, R., Lewandoski, Z., Beyenal, H., 2012