The Rational Design of a Financially Viable Microbial Electrolysis Cell for Domestic Wastewater Treatment

Microbial electrolysis cells (MECs) are yet to achieve commercial viability. Organic removal rates (ORR) and capital costs dictate an MEC’s financial competitiveness against activated sludge treatments. We used numerical methods to investigate the impact of acetate concentration and the distance between opposing anodes’ surfaces (anode apostrophe: interstices width) on MEC cost-performance. Numerical predictions were calibrated against laboratory observations using an evolutionary algorithm. Anode interstices’ width had a non-linear impact on ORR and therefore allowable cost. MECs could be financially competitive if anode interstices widths are carefully controlled (2.5 mm), material costs kept low (£5–10/m2-anode), and wastewater pre-treated, using hydrolysis to consistently achieve influent acetate concentrations >100 mg-COD/l.


INTRODUCTION
Bio-electrochemical systems (BES) promise to revolutionise the wastewater industry by reducing the energy required in wastewater treatment whilst concurrently producing renewable products (Rabaey and Verstraete, 2005). Microbial electrolysis cells (MEC), a type of BES, may pose a viable solution to large-scale energy recovery as they can produce hydrogen gas, a storable form of energy. MEC do not require oxygenated cathode chambers removing the costs associated with aeration. But, possible revenue from hydrogen in MEC over a 20-year period is low relative to the capital costs. Consequently, the largest financial benefit of the technology is in reducing the energy costs of conventional activated sludge (AS) treatments (Aiken et al., 2019). MEC organic removal rates (ORR) must therefore be sufficient to compete with well-established technologies. Our previous financial analyses (Aiken et al., 2019) provided targets for capital costs (c.£600-£1,600/m3.d) and organic removal rates (ORR) (800-1,400g-COD/m3/d) to achieve an MEC that could be financially competitive with AS. Other studies have considered MEC economics but either failed to link costs to ORR (or organic loading rates, OLR) (Sleutels et al., 2012), or calibrated their results against MEC operated at unrealistically warm temperatures (30°C) (Escapa et al., 2012), where microbial activity is a function of temperature (Kolmos, 2012;Lu et al., 2012). This paper seeks to determine whether our targets are achievable, and if so, how. Thereby providing engineers and researchers with a design guide towards a commercially viable MEC for domestic wastewater treatment. Moreover, we seek to do this in a timely fashion. For though MEC have been proposed as a waste treatment technology for nearly 20 years, pilot plant studies are rare and expensive. Exploring performance and design in silico will save both time and money. Previous BES modelling efforts predominantly considered 1D "Monod"-type models for MFC performance (Xia et al., 2018). None have sought to provide a guide towards commercial viability. Such a simulation must be based on realistically available acetate concentrations. Biofilms are dependent upon the mass transport of substrate from the bulk liquid to the anode. Electrogenic species can only, directly, uptake specific substrates (Speers and Reguera, 2012) to fuel their growth and thus electrical generation. Known electrogenic food sources are acetate, formate and hydrogen. Of these, acetate provides the highest current densities when used as a single food source (Speers and Reguera, 2012;Zhao, 2018). As such, acetate concentrations are limiting (Velasquez-Orta et al., 2011;Dhar and Lee, 2014). To provide an empirical baseline for model calibration we observed substrate degradation rates, current response and coulombic efficiencies in a batch-style laboratory set-up fed with realistic acetate concentrations. Our laboratory set-up used controlled conditions, that although idealistic were analogous with the "real-world." For the purpose of simplicity we neglected hydrolysis (and other compounds), but acknowledge that it is limiting. Reactors were fed with a range of acetate concentrations (10-50 mg-COD/l) that reflected the concentration available in domestic wastewater (c.10 mg-COD/ l) (Henze and Comeau, 2008;Huang et al., 2010;Shi, 2011) and a lower range of what may be attainable through hydrolysis (20-100 mg-COD/l) (Ligero et al., 2001).
Varying the space between anodes incurs a trade-off between cost and ORR. For MEC to be viable, we need a decrease in cost and an increase in ORR. Increases in current and hydrogen production rates are secondary (Aiken et al., 2019). Viable large-scale designs often have two anodes facing one another with a "chamber" in-between (Heidrich et al., 2013;Escapa et al., 2015;Cotterill et al., 2017;Baeza et al., 2017), these can vary in width but are usually greater than 10 mm (Supplementary Figure S1A). Alternative porous designs possess anode surfaces that surround a "pore" and similarly oppose one another and are usually less than 1 mm in width (Chong et al., 2019) Figure S1B). We refer to these as the anode "interstices":" the space between anodes' surfaces. As biofilms reside on an anode's surface, packing anode surfaces closer together reduces substrate travel time and thus increases mass transport and subsequently ORR. Though limitations occur at very low interstices' widths (Chong et al., 2019). However, packing anodes closer together also increases anode area requirements. As the majority of capital costs are in the electrodes (Rozendal et al., 2008;Aiken et al., 2019), a costperformance trade-off is incurred. The two key design questions we attempted to answer are therefore: 1) what interstices width maximises ORR whilst minimising cost? 2) at what acetate concentration can MEC be financially competitive with AS?

(Supplementary
We used numerical modelling to obtain MEC costperformance predictions for a range of acetate concentrations and anode interstices widths. We integrated our previous financial viability targets (Aiken et al., 2019) with numerical predictions for ORR and calibrated our models against empirically derived laboratory results. The methods used in this paper were conducted in five phases: "laboratory operation," "biofilm model development," "scheme setup," "calibration" and "extrapolation." To calibrate our numerical model to our laboratory observations we used an evolutionary algorithm, because they are efficient solvers of multi-objective optimisation problems (Weicker and Weicker, 2003). In our final phase, "extrapolation," we consider two scenarios: one in which biomass concentrations were the same as found during calibration and one in which the ultimate biofilm concentration scaled linearly with acetate concentration. ORR outputs were linked with cost-performance targets (Aiken et al., 2019). We varied substrate concentration and anode interstices width, and reported outputs for: ORR, time until 90% substrate removal, anode area requirements and allowable standard costs (/m 2 -anode). We used OpenFOAM as the basis for numerical simulations.

Laboratory Set-Up
Our laboratory set-up used controlled conditions, that although idealistic were analogous with the "real-world." We carefully controlled acetate concentrations (10-50 mg-COD/l), conductivity (770 ± 15μS/cm), temperature (10 ± 0.5°C) and pH (7.50 ± 0.15) to reflect realistic values for an MEC operated with domestic wastewater in a temperate climate. Importantly, acetate concentrations were low to reflect acetate concentrations available in domestic wastewater (c.10 mg-COD/l) (Henze and Comeau, 2008;Shi, 2011;Xin et al., 2018). Additionally, anolyte contained ammonium chloride (20 mg-NH4-equiv/l), a vitamin solution and trace element solution based on wolfe's vitamin solution (ATCC ® MD-VS TM ) and trace element solution (ATCC ® MD-TMS TM ), respectively. We used phosphate buffer (NaH 2 PO 4 , Na 2 HPO 4 ) to control pH and conductivity in both the anolyte and catholyte and a cooled incubator (Samsung, United Kingdom) to control temperature.
We ran eight reactors (Supplementary Figure S2) in "batch" to reduce the impact from convection. Internal compartments for both electrodes were: 40 mm deep, 60 mm wide and 50 mm in height. Reactors each contained a single graphite plate anode (50 mm × 50 mm × 2.5 mm, GPE Scientific UK, machined by Erodex, United Kngdom). Graphite plates were used to create a 'flat' surface. Anodes were placed 35 mm from the outer wall and 2.5 mm from the membrane and sides and thus both sides were exposed to the anolyte. Anolyte volume (V) was 112 ml when displaced by the anode. We glued cost-effective membranes (Rhinohide, Entek, United Kingdom) in between two neoprene frames to provide a watertight seal between chamber compartments. Cost-effectively chosen cathodes (SS316L) were large (A cat 20A an ) so as not to limit current and were placed 2.5 mm from the membrane. We inserted reference electrodes (Ag/AgCl, Cole Parmer, United Kingdom) into salt bridges (1.5% w/v agar (no 1), 1M of KCl) to protect them from the anolyte.

Start-Up
Using potentiostats (Whistonbrook, United Kingdom), we controlled the anode potential of four reactors (−300 mV vs SHE) and applied 600 mV to the other four reactors. We inoculated the reactors with nitrogen sparged settled wastewater taken from a local wastewater treatment works (Chester-le-Street, UK) until we observed the first increase in current. Following this, we fed the reactors with approximately 10 mg-COD/l every 3-4 days. After 4 weeks, 3/8 reactors produced current. After another 4 weeks we changed the feeding window to a 48 h period and 5/8 reactors started up. After another 4 weeks all 8 reactors were producing current. We removed vitamin solution for 2 weeks after we found it to contain lactic acid. Current dropped in all eight of the reactors. Acetate concentrations were then increased to 16.1 ± 3.4 mg-COD/l and a new vitamin solution was administered. Current increased in all 8 reactors and stabilised after 4 weeks.

Operation and Analysis
Following start-up, we fed reactors with sodium acetate every 48 ± 2 h and measured current using the potentiostats (Whistonbrook, United Kingdom). For two separate "feeding runs" during the "plateau phase," we took four anolyte samples (0.5 ml) over a 6 h period (5 min, 1 h, 3 h and 6 h after feeding) to determine acetate degradation rates. We filtered the samples to remove biological material (0.22 μm filter) and flash froze them with liquid nitrogen (−80°C). Upon thawing, we centrifuged samples to remove proteins (10 kDa centrifuge filters, 4°C, confirmed with gel electrophoresis imaging). We measured acetate concentrations (C s ) using Acetate Colorimetric Kits (MAK086-1KT, Sigma Aldrich, United Kingdom) and a 450 nm absorbance spectra (SpectraMax M3, Molecular Devices, United Kingdom). Using SciPy, we determined acetate degradation rates with linear regression (Eq. 1).
q ≈ C s,6h − C s,peak 6h − t peak (1) Once current generation was relatively stable, we observed maximum current (I max ) and current densities (I d, max , Eq. 2, A an 0.005 m 2 ) typically 1-2 h following "feeding" (t peak , Eq. 3). The change in current (current degradation rate) followed an approximately exponential decay curve (Eq. 4). We defined the current degradation rate as the exponent (λ) and computed rates using SciPy. We calculated coulombic efficiencies (Ce peak , Eq. 5) from peak current density and assumed: all substrate (ΔC s ) was removed after 48 h, acetate removal followed linear degradation between the "feeding" time (t feed ) and "maximum current" time (t Imax ), Faraday's constant (96,485C/mol) and the number of electrons supplied from organic mater oxidation (n 8mol-e − /g-COD).

I d, max I max A an
(2) t peak t feed − t Imax (3) Following our "sampling runs," we increased acetate concentrations at a rate of approximately 1.5x the previous concentration (21.8 ± 3.7 mg-COD/l). Actual acetate concentrations varied between reactors and "feeding runs" (Supplementary Table S1). We refreshed media once every 2 weeks and found minor variations in pH (7.50 ± 0.33pH) in both compartments. Three to 4 weeks following an increase in acetate concentration, maximum current stabilised between 'feeding runs'. We repeated our sampling process, increased concentration up to 37.4 ± 9.9 mg-COD/l, waited for current generation to stabilise and sampled again. We used results for acetate degradation rates, maximum current density, current degradation rate and coulombic efficiency from all eight reactors to calibrate our numerical model (Supplementary Table S1). Using these four outputs reduced the degrees of freedom during calibration. There was no statistically significant difference between the systems where anode potential was controlled and voltage applied.

Phase Two: Biofilm Model Development
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 −9 m/s (Leaist and Lyons, 1984). Fluxes (ϕ) were computed using "icoFoam" during "scheme setup" 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 ).
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.
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). 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.

Phase Three: Scheme Setup
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 anodethough 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 2 /s, an estimate for (waste)water at 10°C (Metcalf et al., 2003).

Phase Four: Calibration
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 secondtime 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.
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 fit,Cs f fit + rank f fit , C s i 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).
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 rerunning 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.

Phase Five: Extrapolation
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.

Cost-Performance Analysis
We conducted further analysis to tie our simulated MEC performance with our previous cost-performance targets (Aiken et al., 2019). Our python program computed ORR (Eq. 16) across time and reported ORR at time (τ), the point at which 90% substrate removal was achieved, a target for wastewater removal (European Economic Community Council, 1991).
Our program determined the allowable capital cost per unitflowrate (Cost Q , Eq. 17) using our previous financial model: "target finder" (Aiken et al., 2019). Within "target finder" we removed the costs associated with electrical input and output from the baseline scenario. We equated OLR with ORR and plotted a logarithmic regression between ORR and allowable cost with constants (a,b,c) for 5, 10 and 20 year life-times (n). The program added newly calculated net present costs for electrical input and hydrogen revenue.
To calculate the intermediate hydrogen concentration (C H 2 ) we adapted the equations for current production from organic matter (Logan, 2008) and Faraday's law of electrolysis (Eq.20, n 8mol-e − /g-COD), inputted the change in substrate concentration (ΔC s ) and coulombic efficiency (Ce an 90%), and assumed: cathodic efficiency (Ce cat 50%), molecular weight (M H2 2g/mol) and valency (z 2) of hydrogen. C H2 ΔC s Ce cat Ce an M H2 n z Penultimately, we computed anode area requirements per unit-flowrate (A Q ) from the time until 90% removal (t 90 , Eq. 21), the distance between the centre of the bulk liquid and the biofilm (w int /2), and the thickness of the anode (w an )-equivalent to specific surface area on a "flat" anode.
Finally, we computed standard cost of an MEC reactor (/m 2anode) (Cost A ) as the product of the allowable capital cost per unit-flowrate (Cost Q ) and the inverse of the anode area requirements (A Q , Eq. 22).
We report and compare cost-performance analysis for three acetate concentrations: 10 mg-COD/l, which is reflective of acetate concentrations found in domestic wastewater and two scenarios whose wastewater was assumed to be pre-hydrolysed and achieved influent acetate concentrations of 50-mg-COD/l and 100 mg-COD/l. Outputs are presented as probable outcomes based on the parameters found during calibration. Predictions are described as either the max and minimum or the median plus or minus the interquartile range. Predictions distributions are displayed in figures and are used to measure uncertainty from parameter estimation. Stronger shading relates to more likely outcomes for each scenario and fainter shading relates to less likely outcomes.
Pre-hydrolysing wastewater will allow for financiallycompetitive 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 3 /d) surpassed targets for financial viability (800-1,400g-COD/m 3 /d) (Aiken et al., 2019). And, was more than double the ORR achieved when widths were twice as broad (765 ± 155g-COD/m 3 /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 3 / 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 3 /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 pretreatment, 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 3 /d, 2.5 mm, 100 mg-COD/l) and "linear biofilm development" scenario (1,630 ± 380g-COD/m 3 /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.

MEC Anode Chamber Design
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).
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 2 /m 3 .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.

Standard Unit Cost
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 2 -anode) to break-even are significantly higher in our best scenario (2.5 mm, 100 mg-COD/l, "linear biofilm Frontiers in Chemical Engineering | www.frontiersin.org February 2022 | Volume 3 | Article 796805 development," Figure 3) due to this "double positive" effect. Standard costs to break-even were £16 ± 2/m 2 , £21 ± 2/m 2 and £30 ± 4/m 2 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 , £2.40 ± 1.30/m 2 and £3.60 ± 2.00/m 2 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 2 -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 2 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.

Biofilm Activity
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 2 . Whilst at 1.4 mm current density reduced to 4.7 A/m 2 . 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.

Calibration
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.

Validation
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. FIGURE 5 | Box-plot description for "best fit" parameters following calibration. Largest difference between biofilm and bulk community observed for biomass concentration [1 column].
Frontiers in Chemical Engineering | www.frontiersin.org February 2022 | Volume 3 | Article 796805 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 2 ) better agreed with Chen et al. (2012)'s observations of 7.0 A/m 2 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.

CONCLUSION
Domestic wastewater contains too little acetate (10 mg-COD/l), in its typical form, for MECs to be a viable alternative to AS. ORR are higher when anode interstices are narrower because biofilm activity is greater than bulk community activity. The highest performing MECs are predicted to occur when interstices are sufficiently narrow (2.5 mm). Narrower interstices ( < 2.5 mm) may limit mass transport and thus reduce performance. We estimate that a commercially viable MEC design could be obtained when: standard costs for the new system are reduced to around £5-10/m 2 -anode -wastewater is pre-treated using hydrolysis to obtain a steady supply of acetate ( > 100 mg-COD/l) -and anode interstices' widths (2.5 mm) are carefully controlled. Biofilms must relatively stable over an MEC's operational and quality control must be of a high standard.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.