For our numerical simulation, we adapted the scalar transport equation for diffusion and convection to calculate substrate concentration (C _s ) and included source terms for substrate removal from the biofilm (q _bf ) and the microbial community (q _bc , Eq. 6). Field calculations were conducted within the PIMPLE loop. We set the diffusion coefficient (D _C ) for acetate as 1.2 × 10⁻⁹ m/s (Leaist and Lyons, 1984). Fluxes (ϕ) were computed using “icoFoam” during “scheme setup” d C s / d t + ∇ ⋅ D C ∇ C s ‐ ∇ ⋅ ϕ C s = − q b f − q b c (6)

We added five additional fields within the loop: substrate removal rates from the biofilm (q _bf ) and the bulk community (q _bc ), biomass concentrations of the biofilm (C _bf ) and the bulk community (C _bc ) and volumetric current density (I _v ). Removal rates (q _i , Eq. 8) were computed using Monod kinetics. Monod parameters included terms for the half-velocity constant (Ks _i ) and the maximum specific substrate removal rates (q _i, max, Eq. 7), equivalent to the product of the maximum growth rates (μ _i, max) and the inverse of the yields (Y _i ). q i , max = μ i , max / Y i (7) q i = q i , max C i C s / K s i + C s (8)

Equation 9 linked changes in biomass (dC _i /dt) to removal rates (q _i ), growth yields (Y _i ) and death rates (kd _i ) for both microbial communities. d C i / d t = Y i q i − k d i C i (9)

Adapting Logan (2008) equation gave us volumetric current density (I _v , Eq. 10, F = 96,485C/mol-e⁻, n = 8mol-e⁻/g-COD). We substituted change in concentration for biofilm substrate removal (q _bf ) and assumed that 100% of the electrons released by the biofilm were measured in the cell as current (Ce _bf ). This reduced the degrees of freedom during calibration. In reality, this percentage will be less than 100% and will be dependent upon microbial growth, maintenance and activation losses (Logan et al., 2006), but these contributions are likely to be small (Sleutels et al., 2011). I v = C e b f q b f F / n (10)

To bound the model, we scaled time-step increments in-line with the maximum cell value for substrate removal. No cell within the scheme had greater than 90% removal for any time-step.

In OpenFOAM, we designed our first discretised scheme to reflect a horizontal “cut-out” of our laboratory based anode chamber, using “blockMeshDict” (Supplementary Figures S2, S3A). The scheme consisted of 8 “blocks,” which were situated around the “anode”: represented by boundary walls. We added starting biofilm concentrations to cells adjacent to the anode using “setFieldsDict,” creating a “cross-like” pattern–starting bulk biomass concentration to all other cells–and substrate concentration to all cells. Mesh resolution was 0.5 mm in both the x and y direction. There was only 1 cell in the z-direction.

Following calibration, we created a second discretised scheme. This scheme was made to represent a section of an idealised anode chamber with two ’flat’ anode surfaces on either side of the channel (Supplementary Figure S3B). Alternatively, this could be viewed as an idealised pore channel within a 3D anode - though in reality, anode surfaces are unlikely to be completely flat. The scheme consisted of a single block and the biofilm was situated in a single layer of cells on two opposing sides of the scheme. Again, the bulk community was present in all other cells and starting substrate concentration present in all cells.

Within the “fvSchemes” dictionary: we set time schemes to Euler, gradient schemes to Gauss linear, divergence schemes to Gauss linearUpwind, the Laplacian scheme for substrate concentration to Gauss linear corrected, interpolation schemes to linear and all surface normal gradient terms to orthogonal.

Starting values for fields: biomass and substrate concentrations were controlled within “setFieldsDict.” All other constants and inputs were controlled and inputted in “transportScalarDict.”

Schemes were initialised using “blockMesh.” We set all initial velocities to zero, and pressure across boundary walls to zero. We used “icoFoam,” which uses the “PISO” algorithm to solve the incompressible Navier-Stokes equations for continuity (Eq. 11) and momentum (Eq. 12). The kinematic viscosity (ν) was 1.3 mm²/s, an estimate for (waste)water at 10°C (Metcalf et al., 2003). ∇ ⋅ U = 0 (11) d U / d t + ∇ ⋅ U ⊗ U − ∇ ⋅ ν ∇ U = − ∇ p (12)

Using Python, we created a program to calibrate our model against our empirical observations. The program automatically generated inputs, created our initial scheme, ran our biofilm model, outputted our five fields for each time-step to a csv file and computed model outputs. Computations began from the second-time step to remove simplification errors created during discretisation and initialisation of biofilm concentration values. Inputs included: maximum substrate removal rates (q _bf, max, q _bc, max), death rates (k _d,bf , k _d,bc ), maximum growth rates (μ _bf, max, μ _bc, max), initial biomass concentrations (C _bf , C _bc ), half velocity constants (Ks _bf , Ks _bc ) and yields (Ys _bf , Ys _bc ). Outputs were: substrate degradation rates (Eq. 1), maximum current density (Eq. 2), current degradation rate (Eq. 4) and coulombic efficiency (Eq. 5).

Within this program, we developed an evolutionary algorithm to optimise model inputs (Figure 5) whose resulting outputs were most similar to empirical observations (Supplementary Table S1). The algorithm randomly selected one of the six “feeding run” starting substrate concentrations measured in the laboratory (A-F, Supplementary Table S1). For generation zero–the algorithm generated random values for all inputs for each permutation, within predetermined limits (Figure 5), and ran 100 model permutations in parallel. The algorithm compared model outputs with empirical observations (Supplementary Table S1) using a fitness function (Eq. 13). The fitness included “accuracies,” for substrate degradation rate, maximum current density, coulombic efficiency and current degradation rate. Accuracies were calculated by dividing the output by the mean empirical observation (target) and inversing the result if it was greater than one. We assumed that biomass concentrations were somewhat stable between the start and end of the run due to the “plateau” phase observed in the laboratory. We included two additional parameters into the fitness function to stabilise biomass. The parameters, used biomass concentrations at the second time-step as the target and biomass concentrations in the final time-step as the output. f f i t = ∑ n = 6 x max x n , x = o u t p u t t a r g e t , if o u t p u t t a r g e t < 1 t a r g e t o u t p u t , if o u t p u t t a r g e t > 1 (13)

To reduce overfitting, the function ranked the “normalised scores” in ascending order for each of the six starting substrate concentrations, divided the rank by the total number of iterations (100) and added it to the original fitness score (Eq. 14). f f i t , C s = f f i t + r a n k f f i t , C s / ∑ i (14)

For reporting purposes, mean accuracy was computed as the ratio between model outputs the mean of the four empirical observations: substrate degradation rates, maximum current densities, current degradation rates and coulombic efficiencies (n = 4) (Eq. 15). f a c c u r a c y = ∑ n = 4 x / n (15)

To compute subsequent “generations,” the algorithm selected models with the top 50% of fitness scores. Each selected scheme’s input values was used as a “parent.” The algorithm iteratively selected two “parents” at random and randomly combined their input parameters, using a variable crossover point, to produce a pool of 100 “children.” The algorithm further “mutated” one in six of the “children” by altering one of their 12 inputs at random (±0–2%). The new generation of input permutations were iteratively assigned to a scheme. New substrate concentrations were selected at random and models were run in parallel. This process repeated for 20 generations.

Evolutionary algorithms are stochastic. This means that re-running the calibration will produce marginally different outputs each time. To ensure results were replicable, we ran the model in two batches of 50, compared any differences and took results from the first. To obtain a “global” solution, we selected “best fit” from this plethora of results for the top 10 model permutations for each of the starting concentrations, models also had to have at least 90% similarity between the end and starting biofilm, and outputs had to fall within the ranges observed for maximum substrate degradation rate, current density, coulombic efficiency and current degradation rate. One of the starting concentrations (16.4 mg-COD/l) had no permutations that fell within all four ranges and was excluded from the final result due to incongruency. This process selected a total of 50 input permutations as “best fit.” Viewing the distribution of these “best fit” models gave an estimation of uncertainty from parameter calibration.

We imported each of the “best fit” calibration schemes’ inputs into extrapolation schemes. Our python program varied the width of the channel (2.5–50 mm, 2.5 mm increments) and the starting acetate concentration (10-500 mg-COD/l,10 mg-COD/l increments) before running “blockMesh,” “icoFoam” and the biofilm solver for each scheme. As biomass concentrations are likely a function of acetate concentration, we modelled upper and lower projections for starting biofilm concentrations. Our upper projection assumed that acetate was limiting and that the starting biofilm biomass scaled linearly with starting acetate concentration: “linear biofilm development,” whilst our lower projection assumed that the biofilm concentration remained the same as values found during calibration: “no biofilm development.” Bulk community biomass concentrations remained the same for both scenarios, lest we over-predict substrate degradation rates.

Following model completion, we conducted mesh independence tests to ascertain errors in the discretisation of the scheme. Three cell thicknesses were considered: the original 0.5 mm, as well as two finer resolutions of 0.25 and 0.125 mm. Biofilm thickness was adjusted to be 1 cell thick and the corresponding biofilm concentrations were adjusted to keep biofilm biomass equivalent across comparisons. Discretisation errors were negligible and therefore the models were convergent. Models were also bounded, coulombic efficiency tended towards but was never greater than 100%. However, energy was not 100% conserved by design. Dead microbial cells were not accounted for in our analysis. All “best fit” permutations were stable.

Domestic wastewater contains too little acetate for MECs to be financially competitive with AS (Figure 1). Under the ’simulated domestic wastewater’ scheme (Ac = 10 mg-COD/l), the highest ORR for any scenario (149g-COD/m³/d) was too low to be financially competitive. Higher ORR have been observed in pilot reactors fed with “real” domestic wastewater (130-940g-COD/m³/d) (Cusick et al., 2011; Heidrich et al., 2013; Escapa et al., 2015; Baeza et al., 2017; Cotterill et al., 2017). However, domestic wastewater contains other lipids, carbohydrates and proteins (Huang et al., 2010) that can contribute to ORR but are not directly oxidised by electrogenic species (Speers and Reguera, 2012). Consequently, pilot MECs that did not exhibit hydrogen cycling had low coulombic efficiencies (13–31%) (Heidrich et al., 2013; Baeza et al., 2017; Cotterill et al., 2017). Regardless, empirically observed ORR (Cusick et al., 2011; Heidrich et al., 2013; Escapa et al., 2015; Baeza et al., 2017; Cotterill et al., 2017) are low when compared to AS (450–1,425g-BOD/m³/d) (Stanbury et al., 2017). And lower still than ORR targets for financial competitiveness (800–1,400g-COD/m³/d) (Aiken et al., 2019).

Pre-hydrolysing wastewater will allow for financially-competitive MECs provided anode interstices are sufficiently narrow and biofilms are well-developed. The highest ORR occurred when acetate concentrations were 100 mg-COD/l and interstices were 2.5 mm. Under “linear biofilm development,” ORR (1,630 ± 380g-COD/m³/d) surpassed targets for financial viability (800–1,400g-COD/m³/d) (Aiken et al., 2019). And, was more than double the ORR achieved when widths were twice as broad (765 ± 155g-COD/m³/d). Previously, acetate concentrations of 100mg/L have been achieved by hydrolysing domestic wastewater (Ligero et al., 2001)–though unreliably (20-100 mg-COD/l). ORR are highly sensitive to acetate concentrations (Figure 1). For a more conservative 50 mg-COD/l with an interstices width of 2.5 mm, ORR were 69% lower than achieved at 100 mg-COD/l (508 ± 126g-COD/m³/d, “linear biofilm development,” Figure 1). Accounting for additional substrate removal expected in the “real-world”–this may still be sufficient to achieve viability targets if costs are sufficiently low.

Hydrolysis alone is insufficient for MECs to be viable. Rather, it is the relationship between interstices width and acetate concentration that controls an MECs’ financial competitiveness. MECs can achieve significantly higher ORR when interstices are between 2.5 and 10 mm, under both development scenarios and at all acetate concentrations. Interestingly, this relationship is non-linear. When interstices are greater than 10 mm, interstices width had very little effect on ORR. When anodes were 20–50 mm apart ORR were low (117-185g-COD/m³/d, 100 mg-COD/l). Any successful reactor will therefore need to design anode chambers and treatment systems accordingly. MECs with narrower interstices perform best because organic removal rates (ORR) from biofilm activity are higher than from the bulk community, likely due to higher biofilm biomass concentrations, though we were unable to confirm this due to issues with equifinality. As interstices reduce, more organics are readily available for uptake from the biofilm.

The model failed to compute ORR for the “linear biofilm development” scenario when interstices were less than or equal to 2 mm as the model became unbounded, strongly suggesting that in these narrow chambers there was insufficient acetate to support a sufficiently large biofilm. This was also shown empirically by Chen et al. (2012), who observed MEC current production declining when interstices were less than 2.2 mm and not before. Therefore anodes that have narrower interstices may be limited. Previous pilot reactors have predominantly used “3 dimensional” anode structures such as carbon felt (Heidrich et al., 2013; Escapa et al., 2015; Baeza et al., 2017; Cotterill et al., 2017) and carbon brushes (Cusick et al., 2011). Both materials have high surface area, but also extremely narrow interstices (c.100 μm). The latter may explain why no pilot reactor to date has excelled. If so, novel materials will be needed to achieve suitably high performances.

High quality control during anode manufacturing is needed to ensure all interstices are sufficiently narrow (2.5 mm). However, narrower interstices ( < 10 mm) are likely to result in clogging when fed with “real” domestic wastewater due to the presence of suspended particles (Aelterman et al., 2008). If used for pre-treatment, up-flow hydrolysis reactors can retain some, but likely not all, suspended solids (SS). Further pre-treatment, such as screening, will be required to remove and dispose of SS, increasing cost. Wider interstices (10 mm) would mitigate this problem, but would significantly reduce ORR (78% lower). To facilitate sufficient mass transport in narrow interstices large MEC may require a modular design, increasing cost. Alternatively, convection may be increased, though an increase in turbulence could reduce biofilm size (Fink et al., 2016).

If biofilm development is poor then ORR will not reach targets. Poor biofilm development had the highest impact on ORR. The difference in ORR between the “no biofilm development” scenario (473 ± 141g-COD/m³/d, 2.5 mm, 100 mg-COD/l) and “linear biofilm development” scenario (1,630 ± 380g-COD/m³/d, 2.5 mm, 100 mg-COD/l) is highly significant. As biofilm development likely lies between our two projections, further empirical work will be needed to confirm the results from these, though we hope that these projections provide researchers and technologists with a floor and ceiling for MEC performance, in temperate conditions, with respect to acetate concentrations and interstices’ widths. Moreover, biofilms typically undergo growth and death cycles and further work will also be needed to investigate the risk of biofilm failure.

Faster reactors are smaller reactors. The colonisable anode surface area needed to support sufficient biofilm development to treat a given flowrate is a function of the interstices width and the size of the reactor (Figure 2). Under “linear biofilm development,” when widths are broad ( > 12 mm) and concentrations are high, time taken to reach 90% substrate removal is longer. However, when widths are narrow ( < 12 mm) and concentrations are high, time taken to reach 90% substrate removal is shorter. Narrower interstices ( < 12 mm) therefore reduce both an MEC’s retention time, size and anode surface area when acetate concentrations are high (100 mg-COD/l).

The fastest, and thus the smallest reactor (100 mg-COD/l, 2.5 mm, “linear biofilm development”) was predicted to achieve 90% removal of acetate in 1h18 ± 0h23 and will require 21.7 ± 6.4m²/m³.d to support “sufficient” biofilm biomass. If interstices are narrow but pre-hydrolysis is poor (50 mg-COD/l, 2.5mm, “linear biofilm development”) 90% removal will be marginally slower (2h06 ± 0h30) and larger anodes are needed to support a more diffuse biofilm (35.1 ± 8.4m²/m³.d). If anode quality control is poor (100 mg-COD/l, 5 mm, “linear biofilm development”) then removal will be significantly longer (2h49 ± 0h35)–though area requirements are only marginally higher (23.4 ± 4.9 m²/m³.d). Furthermore, if interstices are wide, to prevent blockages (100 mg-COD/l, 10 mm, “linear biofilm development”), then retention times will be significantly larger (6h05 ± 1h21), and increases in anode area are needed to treat the same flowrate (45.8 ± 8.7m²/m³.d). Most importantly, biofilms must be well developed and relatively stable. Poor biofilm development (100 mg-COD/l, 2.5 mm, “no biofilm development”) will require significant increases in both removal time (4h35 ± 1h16) and anode area (76.4 ± 21.0m²/m³.d).

If interstices are narrow ( < 10 mm) and acetate concentrations high (100 mg-COD/l), MECs could have significantly smaller HRT than typical for AS (5–14 h) (Stott, 2003). Therefore, an MEC would be retrofittable into aeration chambers, provided cathodes are small. However, this presupposes that most of the COD in the influent is converted into acetate or other quickly digested VFAs. Failing this, MECs will require longer retention times or cheaper post-treatments to “polish” effluent to standard. Moreover, removal times will likely vary with substrate concentrations and biofilms undergo growth and death cycles (Judd, 2013) and therefore larger anode areas are needed. Considering biofilm stability and applying a safety factor around 2–5x our ideal scenario, a rough conservative measure for anode area per flowrate could be 50–100 m²/m³.d, depending on the material, cost, prior testing and risk tolerance. Material surface area should only count towards area requirements if the surface is accessible and capable of supporting an electroactive biofilm. Materials with very high specific surface areas, such as activated carbon cannot have the majority of their surface counted, as the high surface area is predominantly in the form of a multitude of pores with small diameters ( < 50 μ m). The most suitable materials will be flat carbon surfaces that can be packed tightly together and porous carbon materials with pores greater than 2 mm. We challenge researchers and technologists to acquire or create new anodes that fit this requirement.

Reactors with densely packed anodes benefit from higher ORR and thus higher allowed capital costs (Aiken et al., 2019). Moreover, sufficient anode area requirements reduce when interstices are narrow and concentrations are high (2.5 mm, 100 mg-COD/l, “linear biofilm development” Figure 2). Consequently, allowable standard costs for MEC reactors (standardised/m²-anode) to break-even are significantly higher in our best scenario (2.5 mm, 100 mg-COD/l, “linear biofilm development,” Figure 3) due to this “double positive” effect. Standard costs to break-even were £16 ± 2/m², £21 ± 2/m² and £30 ± 4/m² for 5, 10 and 20 year time-frames, respectively. However, if biofilm development is poor (2.5 mm, 100 mg-COD/l, “no biofilm development”) standard costs are £1.60 ± 1.04/m², £2.40 ± 1.30/m² and £3.60 ± 2.00/m² for 5,10 and 20 years, respectively.

Biofilm growth was the largest factor affecting standard costs (Figure 3). If biofilm development is close to the “linear biofilm development” scenario then allowable standard costs appear sufficiently high that future MEC technologists should be able to work with materials whose prices are within these bounds. However, if biofilm development is closer to the “no biofilm development” scenario It is unlikely that sufficiently cheap materials can be sourced: the cheapest material in a recent MEC pilot was the membrane (£1/m²-membrane, Rhinohide, Entek, United Kingdom) (Aiken et al., 2019).

In reality it is most likely that mean performance will be somewhere between the “linear biofilm development” and “no biofilm development” scenarios, and removal will not consistently meet standards (European Economic Community Council, 1991). Reducing material costs to around £5–10/m² would allow for some tolerance in development, whilst still retaining the possibility of being cost-effective in the long-term. Though, this will not protect against completely poor performance. Researchers and technologists will need to base any costs on their own reactor’s performance.

ORR is higher with narrower interstices’ because biofilm activity is greater than bulk community activity (Figure 5). As interstices’ become narrower, the biofilm plays a larger role in substrate degradation as shown by the increase in coulombic efficiency (Figure 4). There is insufficient substrate available to replace that which is consumed by the biofilm. Maximum current density remains the same at narrower interstices but substrate and current both degrade faster (Figure 4). This was true of both scenarios, regardless of parameter uncertainty.

Our model assumes that biomass concentrations remain the same regardless of width. As current density is a function of the biofilm, we believe this serves as a proxy. Chen et al. (2012) observed that in an MFC fed with 20 mM of acetate (1,260 mg-COD/l) current densities declined when interstices’ widths were reduced below 2.2 mm and not before. At interstices widths of 2.2, 2.9 and 3.2 mm current densities were around 7.0A/m ². Whilst at 1.4 mm current density reduced to 4.7 A/m². Chong et al. (2019) conducted a review of 71 papers and found that the greatest reduction in planar current densities occurred when pore sizes were <0.5 mm due to pH inhibition. This was further speculated as an important limiting factor by Logan et al. (2006). However, only 5 out of the 71 papers had anode interstices with widths >0.5 mm and broader chambers (>1 mm) were not considered.

We did not consider the impact of a reduced HRT on biofilm development. If all substrate is consumed in a short period of time reactors could be fed more frequently, and could hypothetically support a larger biofilm as more substrate would enter the biofilm. However, this would only further highlight our findings on the importance of the relationship between anode interstices’ width and cost-performance. We encourage experimentalists to investigate further.

In multi-objective problems, there exist a set of optimal solutions that trade off between different objectives, known as the Pareto Front (Weicker and Weicker, 2003). Obtaining higher accuracies beyond this front is impossible as any improvement in one objective leads to a deterioration in one or more of the other objectives. The two best fitting schemes had a “mean accuracy” of 96.4% with starting acetate concentrations of 21.2 mg-COD/l and 39.4 mg-COD/l. The best fitting schemes for the remaining concentrations: 14.7 mg-COD/l, 16.4 mg-COD/l, 22.4 mg-COD/l and 35.4 mg-COD/l, had mean accuracies of 87.6%, 88.4% 94.6 and 86.9%, respectively. As the algorithm was not designed to map the whole parameter space, it is possible that better fitting parameters exist for these concentrations but that the fitness function was unable to select for them. However, as accuracy was measured against a particular concentration, a higher accuracy does not necessarily mean that the parameters of that run fit well with the other concentrations.

Moreover, empirical observations are themselves a model. Measurements were made using a combination of chemical and electronic signals that present information about an underlying reality but do not exactly map it. Uncertainties in empirical measurements mean that the lab results themselves cannot be entirely relied upon. Indeed, one of the starting concentrations (16.4 mg-COD/l) failed to find any combination of parameters that had outputs that were within all four of the observed ranges. This is likely because the feeding run at 16.4 mg-COD/l had a lower maximum current density and coulombic efficiency than the feeding run at 14.7 mg-COD/l, despite containing a higher starting acetate concentration and having a larger substrate degradation rate.

We present our results to demonstrate uncertainties resulting from parameter estimation of mean observations for microbial electrolysis cells with well-developed biofilms. Uncertainty increased as concentration increased. In the laboratory we fed our reactors with low acetate concentrations (10-50 mg-COD/l). As a result, the model was not constrained above these concentrations and as the concentration of acetate became greater than 50 mg-COD/l the predictions from the model were more uncertain due to uncertainties in the parameters (Figure 5). By uniformly randomising the starting concentrations for each iteration within each of the generations, the algorithm was only softly constrained by any one measurement. Permutations were assessed, each generation, based on their selective fit against any one of the concentrations and in the subsequent generation the parameters were rearranged and assessed against any other. If one concentration, or a group of concentrations, fit better than another, then the permutations that had poorer fit would be selected against, thus reducing overfitting. However, due to the stochastic nature of the algorithm, there still lies a small possibility, that a high accuracy measurement fits one of the starting concentrations very well but not the others. Therefore, observing the extrapolation results as probabilistic and viewing the median gives a more accurate picture of the “global” solution. Furthermore, rerunning the algorithm many times further reduces this possibility and the impact that random starting permutations, and combinations have on the final output. By looking at the top results for each starting concentration, rather than the most accurate solutions, we consider both the global solution and the uncertainties in parameter estimation.

Parameter estimation was further confounded by issues with equifinality in some of the parameters. A notable example is that of maximum substrate degradation rate and biomass. The substrate degradation, which we calibrated for took the product of these parameters (Eq. 1) and neither was constrained by observation. Extrapolating to higher concentrations beyond 100 mg-COD/l showed a wider range of predicted substrate degradations and current densities. Yet, despite this solutions converged when interstices widths were narrower. Moreover the non-linear relationship between width and ORR held for all permutations regardless of whether the concentration of biomass increased or remained stable with changes in acetate concentration. This is because in order for the mechanics of the system to fit with observed current densities and substrate degradation rates microbial activity in the biofilm must have been higher than in the bulk liquid.

None of our observations particularly favoured either the “no growth” or “linear growth” scenarios when plotted against acetate concentration (Figure 6). Moreover, very few studies investigate electrogenic biofilm activity at low temperatures despite activity being a function of temperature, and those that do consider different substrates at much higher COD concentrations (Larrosa-Guerrero et al., 2010; Lu et al., 2012). In all likelihood performance lies somewhere between these two scenarios and we present these as a floor and a ceiling for MEC performance.

Our model predicted substrate concentration well, between 10 and 50 mg-COD/l. Both scenarios predicted similar substrate degradation rates at these concentrations. Beyond 50 mg-COD/l the scenarios began to diverge. At 100 mg-COD/l the “linear biofilm development” predicted higher substrate degradation rates (6.3 ± 2.7g-COD/l.h vs 5.4 ± 2.9g-COD/l.h). Higher acetate concentrations showed that the “linear biofilm development” scenario continued to achieve higher substrate degradation rates. This is unlikely in reality as higher concentrations can have an inhibiting effect (Mateo et al., 2019). We therefore believe that the “linear growth” scenario over predicts substrate degradation rate when concentrations are high (>100 mg-COD/l) and interstices’ wide (>10 mm). Moreover, Zhao (2018) found that when reactors were acclimated with high acetate concentrations Geobacteraceae dominated (91% of anodic microbial community, qPCR). Thus, bulk community biomass concentrations may diminish as electrogenic organisms outcompete them for substrate.

Maximum current density is a function of the biofilm, whilst substrate degradation rate is a function of both the biofilm and bulk community. Therefore ORR at narrower interstices’ when the biofilm dominates is most affected by our prediction of maximum current density. Divergences in maximum current density between the two scenarios are stark (Figure 6). Observations for maximum current density had a linear relationship with acetate concentration, as was expected from observations. However, our observations appear to lie between the “no growth” and “linear growth” scenarios. Though there is far too much variability in measurements and too little difference between the scenarios at low concentrations to give credit one way or the other using these measurements. At 2.5 mm, 1,260 mg-COD/l, our “linear biofilm development” scenario (9.5 ± 12.9 A/m²) better agreed with Chen et al. (2012)’s observations of 7.0 A/m² than our ’no biofilm’ development scenario (0.2 ± 0.2 A/m2). However, uncertainties in our predictions were high relative to the median and this experiment was conducted at a higher temperature. Further work with low temperatures (10°C) and narrow anode interstices’ (c.2.5 mm) is needed.

Coulombic efficiency was not well mapped by changes in acetate concentration (Figure 6). In the lab we observed a large range of efficiencies at relatively small concentrations, which may mean that that the link between these phenomenon is either negligible or indirect. Current degradation rate’s relationship with acetate concentration appeared to be best mapped by the “no growth” scenario (Figure 6). However, the relationship between current degradation rate and acetate concentration is still unclear and is not well understood. Furthermore, the “no growth” scenario became somewhat unstable at acetate concentrations greater than 50 mg-COD/l. Whilst, the “linear growth scenario” remained stable and similarly mapped the centre of the data points we observed (Figure 6). Further study is needed before any conclusions on current degradation rate can be made.