Vincent Feynerol
email: vincent.feynerol@univ-lorraine.fr

Ranine El Hage

Mariela Brites Helú

Vanessa Fierro

Alain Celzard

Liang Liu

Mathieu Etienne
email: mathieu.etienne@univ-lorraine.fr

Comparative kinetic analysis of redox flow battery electrolytes: from micro-fibers to macro-felts

Keywords:

In the context of new aqueous organic electrolytes development for redox flow batteries, two different approaches of charge transfer kinetic constant determination on carbon-based felts are compared. The first one consists in direct electrochemical measurement on centimetric pieces of felts, referred to as macro-felts, and the second is based on electrochemical measurement on carbon fibers extracted from felts and conditioned as microelectrodes, referred to as micro-fibers. Quantitative cyclic voltammetry (CV) involving experimental data fitting or simulation and electrochemical impedance spectroscopy (EIS) are carried out on the two supports, with two typical aqueous organic electrolytes, and determined values for charge transfer kinetic constants are compared. k 0 values analyzed from micro-fibers are systematically lower than that analyzed from macro-felts where the differences may reach ca.

two orders of magnitude in some cases. CV and EIS analysis on micro-fibers show very good consistency. With the simplicity of physical analysis, micro-fibers are considered as a reliable and powerful tool for kinetic determination, especially when confronted to fast kinetics.

Introduction

Redox flow batteries are a promising energy storage technology that may play an important role in the integration of intermittent renewable energies in grid networks. Typical electrodes used in redox flow batteries are porous materials that offer high specific surface area. Carbonbased felts are commonly used as electrode material in redox flow batteries due to their high specific surface area, high conductivity, resistance to corrosion and wide operation potential range in addition to their relatively low cost [START_REF] Liu | Progress on the electrode materials towards vanadium flow batteries (VFBs) with improved power density[END_REF].

Although carbon felts allow easy permeation of electrolyte, yielding high current density in practical applications, their intrinsic electrochemical behavior is not easy to analyze. The apparent electrode performance depends not only on the properties of the electrode material but also on their geometry which is often not well defined for porous felts. Most electroanalytical methods are based on regular geometry (planar, cylindrical, spherical, microdisk, etc.), and model equations are usually developed for planar electrodes. In the case of porous electrodes, the mass transfer within their structure shall be considered, which greatly complicates the mathematical analysis. The high surface area also makes the differentiation of apparent and intrinsic charge transfer kinetics difficult, especially when the latter could also be very fast in most energy storage systems. On the other hand, for guiding the material development, it is always of high interest to assess the intrinsic charge transfer kinetics of porous carbon felts.

Cyclic voltammetry (CV) and electrochemical impedance spectroscopy (EIS) are the most common electrochemical methods to determine charge transfer kinetics. The CV analysis for planar electrodes is highly matured [START_REF] Compton | Understanding Voltammetry[END_REF] and one may derive kinetic information by routinely fitting of experimental data with appropriate models using commercial software. For analyzing porous electrodes, the geometry needs to be properly approximated. One approach is considering arrays of electrode units, with the geometry chosen according to the type of porous material studied [START_REF] Barnes | Voltammetry at porous electrodes: A theoretical study[END_REF], [START_REF] Smith | The electrochemical characterisation of graphite felts[END_REF]. Very recently, Tichter and co-authors improved this formalism by representing porous electrodes as randomly distributed arrays of microelectrodes for wider variety of geometries [START_REF] Tichter | Theory of cyclic voltammetry in random arrays of cylindrical microelectrodes applied to carbon felt electrodes for vanadium redox flow batteries[END_REF], [START_REF] Tichter | Universal Algorithm for Simulating and Evaluating Cyclic Voltammetry at Macroporous Electrodes by Considering Random Arrays of Microelectrodes[END_REF]; they developed a software, Polarographica, for simulating CV of porous electrodes.

Considerable work has also been made to analyze kinetics from EIS of porous electrodes.

DeLevie is often referred as pioneer of the field in the 1970s [START_REF] De Levie | On porous electrodes in electrolyte solutions. I. Capacitance effects[END_REF], [START_REF] De Levie | On porous electrodes in electrolyte solutions-IV[END_REF] with the introduction of transmission line analogy. Following this work, Candy and co-authors studied gold powders' impedance in KOH [START_REF] Candy | The characterization of porous electrodes by impedance measurements[END_REF], [START_REF] Candy | The pore texture of raneynickel determined by impedance measurements[END_REF] and obtained good fit from numerical calculations considering their catalyst as an array of identical cylindrical pores. In the 1990s, Lasia established a similar numerical model including a potential gradient within a porous electrode and later added a concentration gradient [START_REF] Lasia | Impedance of porous electrodes[END_REF]- [START_REF] Lasia | Modeling of Impedance of Porous Electrodes[END_REF]. Paasch and co-authors developed a "macrohomogeneous material" model [START_REF] Paasch | Studies on thermally structurized polyacrylonitrile by the electrochemical impedance method[END_REF], [START_REF] Paasch | Theory of the electrochemical impedance of macrohomogeneous porous electrodes[END_REF], and similar approaches were further complexified in cases of problems related to Lithium rechargeable batteries [START_REF] Meyers | The Impedance Response of a Porous Electrode Composed of Intercalation Particles[END_REF], [START_REF] Doyle | Computer Simulations of the Impedance Response of Lithium Rechargeable Batteries[END_REF].

As the challenge of directly analyzing porous felts from both CV and EIS mainly arises from the undefined geometry, an alternative way is to extract the fibers from the felt and study them as microelectrodes. As the diameter of fiber is in the range of a few microns while the length of fiber is usually in millimeter range, the contribution of the tip can be neglected and the electrode can be regarded as a 1D cylindrical system which eases the mathematical treatment. Kovach et al. used such electrodes with various electrochemical systems and analyzed in particular their results by fitting experimental CV [START_REF] Kovach | Electrochemistry at partially blocked carbon-fiber microcylinder electrodes[END_REF]. Bourke et al. [START_REF] Bourke | Effect of Cathodic and Anodic Treatments of Carbon on the Electrode Kinetics of V IV /V V Oxidation-Reduction[END_REF] and Miller et al. [START_REF] Miller | Kinetic Study of Electrochemical Treatment of Carbon Fiber Microelectrodes Leading to In Situ Enhancement of Vanadium Flow Battery Efficiency[END_REF] carried out CV and EIS on carbon fibers to measure the effect of anodic and cathodic treatments on charge transfer kinetics with typical vanadium electrolytes. Recently, Landon-Lane et al. carried out EIS on carbon fiber microelectrodes with ferrocyanideferricyanide and proposed a simplified approach to fit data for cylindrical geometry [START_REF] Landon-Lane | EIS at carbon fiber cylindrical microelectrodes[END_REF]. The author further developed the use of microelectrodes by fitting voltammograms obtained with both ferrocyanide-ferricyanide and vanadium electrolytes [START_REF] Landon-Lane | Investigating Electron Transfer Kinetics in Redox Flow Batteries[END_REF], minimizing the difference in peak potential separation between experimental and simulated data, conforming to methods suggested by other authors [START_REF] Aoki | Theory of irreversible cyclic voltammograms at microcylinder electrodes[END_REF], [START_REF] Neudeck | The determination of diffusion coefficients and rate constants from the dependence of the peak separation and peak current on the scan rate of cyclic voltammograms at micro-cylindrical electrodes[END_REF].

In this work, we carried out comprehensive studies of intrinsic charge transfer kinetics of two typical electrolytes promoted in the development of aqueous organic electrolyte redox flow batteries on activated and non-activated carbon felts: ammonium ferrocyanide and a derivative of TEMPO. The analysis was achieved by fitting CV and EIS data on both macrofelts and extracted micro-fibers and the results were compared to identify the best experimental and modelling approach to quantify charge transfer kinetics of redox flow electrolyte on such porous material.

Materials and Methods

Chemical products

Two aqueous organic electrolytes are presented in this study: ammonium hexacyanoferrate, received as pure ammonium ferrocyanide (NH 4 ) 4 Fe(CN) 6 , and a derivative of TEMPO (one of the cycle's hydrogen is substituted with an organic group), further referred to as "TEMPObased molecule", both provided by CIC energiGUNE. Solutions were prepared by dissolving the right amount of salt in its corresponding supporting electrolyte: 1 M NaCl for TEMPObased molecule, prepared from pure NaCl (Prolabo, France), and 1 M NH 4 Cl for ammonium hexacyanoferrate, prepared from pure 99.5% NH 4 Cl (Fisher Scientific, USA). The concentrations of the studied solutions were 10 -4 M for "macro-felts" and 10 -2 M for microelectrodes.

Experimental apparatus and protocols for macro-felts

All electrochemical experiments were performed with a VSP-300 potentiostat (Bio-Logic, France) interfaced by EC-Lab® software (Bio-Logic, France).

For "macro-felts": a custom-made PVC cell to perform cyclic voltammetry and electrochemical impedance spectroscopy on pieces of graphite felts was designed and built (see Figure S1). Graphite felts were used as electrodes (SGL 4.6 GFD, thickness of 4.6 mm): two large cylindrical pieces (Ø 5 cm) served as counter electrodes (CE) while one smaller cylindrical piece (Ø 1 cm) was the working electrode (WE). The felts were cut with punches for better reproducibility. Both commercially available non-activated (n.a.) and activated versions of this felt were tested. The electrodes were positioned inside the cell using 3Dprinted frames so that the WE was centered between the two CE and set so their respective centers were aligned with one another on a horizontal line. A 3D-representation in Figure S1 shows the positioning of the electrodes. Graphite felt electrodes were connected to the leads of the potentiostat by piercing through each of them with a composite rectangular carbon rod pultruded (DPP™, R&G Faserverbundwerkstoffe GmbH) attached to banana connectors with copper wire and silver-epoxy resin (silver conductive epoxy adhesive, MG Chemicals). The reference electrode was chosen to be an Ag/AgCl quasi-reference electrode (QRE) in order to mitigate interference with the measurements. QRE was made by coating a Ø 0.25 mm silver wire (99.9%, chemPUR Feinchemikalien und Forschungbedarf GmbH) with AgCl by applying 1 V between the wire and a large piece of graphite felt as the counter electrode in 1 M HCl for one minute. This wire was then inserted in the cell via one of the frames to avoid contact with any of the other electrodes. The cell was closed with a 3D-printed lid. Before each experiment, felt cylinders were hand pressed in the studied electrolytic solution in the manner of a sponge in order to expel the air contained in the pores and to fill them with electroactive solution, as non-activated materials have poor water wettability.

EIS experiments for graphite felts were carried out applying the following protocol: first, the formal potential of studied redox couple was estimated from cyclic voltammetry by taking the average of oxidative and reductive peak potentials. Then, staircase potential electrochemical impedance spectroscopy (SPEIS) was performed with EC-Lab® software in a range of potentials centered on formal potential and including it: at each potential step, potential is applied for enough time to reach pseudo-equilibrium (typically ten minutes), then oscillates around its pseudo-equilibrium value while measuring the response current at different frequencies. For ferrocyanide, 50 mV potential steps and five minutes of buffer time were used. For TEMPO-based compound, 20 mV potential steps were applied using a buffer time of two minutes. In both cases, oscillations amplitude was set at 10 mV. Screening different potentials allowed us to verify that a minimum impedance indeed occurred at formal potential for our two electrolytes, therefore validating the assumption of a Nernstian system.

Performing the experiments around formal potential rather than equilibrium potential (or open-circuit potential), as often found in literature, was motivated by the fact that the TEMPO-based molecule received was only present in its reduced form and was extremely stable: it was thus impossible to use i-E equations for analysis as the oxidant surface concentration was unknown, and could not be approximated to 0 since the system's kinetics were very fast compared to mass transfer. Frequencies screened ranged from 7 MHz to 50 mHz, with six frequencies screened per decade, and for a total impedance measurement duration of 2.5 minutes.

The analysis of EIS spectra was done by performing Complex Non-linear Least Squares, that is to say by fitting the experimental Nyquist diagram with a model. The models used are detailed in a subsequent part. The optimization was done by programming a small custommade routine with MATLAB® software based on the use of the optimization function "fmincon". Constraint tolerance, function tolerance and optimality tolerance were all set to 10 -8 . The maximum number of function evaluations was set to 200 and algorithm used was "sqp". Before optimization, the parameters were weighted with their approximate order of magnitude to ease convergence. The relevant range of frequencies usable for analysis was typically below 5 kHz for graphite felts.

BET surface of graphite felts was measured using 3 g of graphite felt. The reproduction of the experiment showed a variation around 10% of estimated value.

Experimental apparatus and protocols for microelectrodes

For microelectrodes: the cell used to perform EIS and CV experiments was made out of a simple flat-bottomed polystyrene sampling tube (Ø 3.5 cm) closed with a 3D-printed lid. The lid allowed to mount and fix the position of an Ag/AgCl QRE (prepared similarly to the one used with graphite felts), a counter electrode made from a DPP™ composite rectangular carbon rod pultruded and a custom-made 3D-printed support holding the carbon fiber microelectrode studied. The 3D-printed support for fiber microelectrodes consists in two parts: part A, longer, serves to allocate a carbon fiber and a copper wire (Ø 0.2 mm), and part B, shorter, seals the electrical contact fiber-copper and completes the support. The groove in part A is filled with silver-epoxy resin (Silver conductive epoxy adhesive, MG Chemicals).

From one side, copper wire is inserted into the guiding hole until it reaches the Ag-epoxy. On the other side, carbon fiber is position on the groove ensuring the good contact with the copper wire. The piece is cured at 65°C for 2 h. Then, the exposed electrical contact is sealed by gluing part B and part A using UV-curing resin. The microelectrodes prepared as described show excellent chemical and mechanical stability over time in the tested electrolytes.

Photograph as well as an exploded view are available in Figure S2 for better visualization.

Before each experiment, microelectrodes were cleaned by immersing the tip of the support in a 70% ethanol solution for ten minutes in order to remove possible presence of dust or other impurities on the surface of the fibers.

Carbon fibers' mean radius was estimated after observation of several carbon fibers with scanning electron microscopy (JSM-IT500HR, JEOL). They were found to range typically between 8 and 10 µm, thus an average diameter of 9 µm was chosen, therefore giving a mean radius of 4.5 µm. Carbon fibers' lengths were estimated by optical inspection under the microscope, with an estimated error of ± 0.25 mm.

The procedure applied for EIS experiments with carbon fiber microelectrodes was the same as the one with graphite felts, except for ammonium hexacyanoferrate on non-activated fiber: in this particular case, the kinetics were found very slow and an impedance minimum could no longer be found at formal potential, instead appearing at more oxidative potential. In this case only, a wider potential range was scanned, and the analysis was focused around the potential giving minimum impedance values. This choice was motivated, first by the impedance spectra shape at these potentials being more consistent with the one theoretically obtainable with the model used, and second by failure to obtain good fit with said model at equilibrium potential.

In addition, due to the very small surface of microelectrodes, frequency could be analyzed up to very high frequencies, typically up to 200 kHz.

All CV and EIS experiments were triplicated for reproducibility.

Mathematical modelling 3.1 "Macro-felts" porous electrodes modelling

Cyclic voltammetry model for porous electrodes

In order to estimate kinetic constants based on cyclic voltammograms measured on porous macro-felts, the model developed by Tichter and co-authors was used. This model, available under the form of a free software called Polarographica, is based on the representation of porous materials as randomly spaced arrays of microelectrodes. The spacing of microelectrodes is handled with a statistical distribution that can be divided into several intervals of inter-electrode distances. Cyclic voltammetries are then calculated for each of these intervals, each time based on corresponding average distance, and their sum is weighted according to the statistical distribution to give the resulting voltammogram. Among the different geometries available, the array of cylindrical electrodes is the one that should best represent carbon and graphite felts. Although currently unable to automatically fit curves, it is possible to import experimental data and to manually change input parameters so for a simulated curve to approach the experimental one.

Using version 2.1.0 of this software with "statistically weighted cylindrical, cylindrical finite diffusion model", CV curves were simulated with typical parameters determined by the authors and obtained with similar felts [START_REF] Tichter | Universal Algorithm for Simulating and Evaluating Cyclic Voltammetry at Macroporous Electrodes by Considering Random Arrays of Microelectrodes[END_REF]. Hence, the distribution of inter-electrode distances was not changed from "default parameters", and intervals of distribution as well as the number of cylinders per mm 2 were set to the same values as the ones used by the authors with the same type of felt (respectively 722 cylinders per mm 2 and 10 intervals). Then, the total surface was adjusted so for the peak magnitude to match the experimental one and the kinetic constant was adjusted in order to find the minimum value allowing the peak-to-peak distance to match the experimental one.

Macrohomogeneous model for graphite felts EIS

For analysis of the impedance spectra with graphite felts, the macrohomogeneous model developed by Paasch, Micka and Gersdorf [START_REF] Paasch | Theory of the electrochemical impedance of macrohomogeneous porous electrodes[END_REF] was used, in addition to some minor modifications. This model assumes the porous electrode as a biphasic material, composed of a liquid phase (the electrolyte inside the pores) and a solid phase (the actual electrode). The porous electrode is in contact on one side with the bulk solution and on the other side by a metal connector. This model is one-dimensional: physical variables are assumed constant over one "slice" of the porous element. The concentration of active species is assumed constant through the whole electrode, but the model allows the addition of a "hindrance function" in order to represent for instance the influence of diffusion phenomena. Finally, double-layer capacitance per length of porous material is also assumed constant. This model considers the separation of the total electric current circulating through the electrode to the part circulating in solid phase and to the part circulating in the liquid phase. Their respective value changes with the axial position in the porous electrode, but their sum stays constant by virtue of charge continuity. A gradient of potential with position in solid and liquid phase therefore arises. This model was originally derived relative to equilibrium potential, but it can be shown that the equation they established is basically identical when formal potential is used instead of equilibrium potential (development in supplementary material). The main difference is situated in the expression of parameter . The complex impedance equation is given in equation ( 1), the parameters are further explained in the following equations and symbols are given in appendix:

(1) (2) (3) (4) (5) (6) (7)
R ohm is the ohmic resistance between the working electrode and the reference electrode.

n is the number of electrons involved in the electrochemical process (equal to 1 for both salts), F is the Faraday's constant, R g is the ideal gas constant, T the room temperature (298.15 K). d is the electrode thickness, equal to 4.6 mm for all cases studied her. ρ 1 and ρ 2 are the resistivities of respectively the solid phase (the actual electrode material) and the liquid phase (the electrolyte) within the porous material. ρ 2 is modified from the bulk electrolyte resistivity ρ 2 * due to the porosity ε and tortuosity f p of the material according to equation ( 10):

(10)
k 0 is the intrinsic charge transfer constant while is the specific surface area of the porous electrode, that is to say its total reactive surface relative to its geometric volume. A is the geometric area of the porous electrode cross section.

Parameter β was modified to include a constant phase element (CPE) instead of a capacitor to represent the double-layer of the porous material. This can be done simply by adding an exponent as shown in equation [START_REF] Compton | Understanding Voltammetry[END_REF]. The capacitance C in parameter k and K then becomes a pseudo-capacitance.

The mass transfer "hindrance function" was chosen to represent finite diffusion with reflective boundary. The finite diffusion reflective boundary, as reported by MacDonald [START_REF] Macdonald | Impedance Spectroscopy -Emphasizing solid materials and systems[END_REF],

corresponds to solving the diffusion equation with a specific boundary condition. The use of this boundary rather than a transmissive boundary was mainly motivated by the shape of our impedance spectra that tended towards infinite values at low frequencies for the imaginary part of the impedance. In addition, a modification to the finite diffusion hindrance function was performed in order to correctly represent our data at low frequencies: a dimensionless exponent ϕ was added to represent "non-uniform diffusion" [START_REF] Macdonald | Impedance Spectroscopy -Emphasizing solid materials and systems[END_REF]. This modification allows the slope in the Nyquist diagram to tend towards a finite value instead of infinity, which was observed in particular for ammonium hexacyanoferrate. A similar modification was performed by Sun et al. [START_REF] Sun | Resolving Losses at the Negative Electrode in All-Vanadium Redox Flow Batteries Using Electrochemical Impedance Spectroscopy[END_REF], although the parameter introduced was a factor rather than an exponent. The values of fixed physical parameters taken as input in the model are given in Table 1, parameters fitted were: R ohm , C, γ, k 0 , l p , S c , f p , ϕ.

The CPE exponent was not always needed for the convergence to occur, and it was fixed to 1 unless convergence was impossible without it. Failed convergence would typically manifest in the routine not being able to find a value for the kinetic constant, and constantly searching for a higher value, stopping at the upper boundary. Increasing the value of this upper boundary would result in the model reaching a higher value for k 0 while leaving the other parameters virtually unchanged. In these cases, the value found for ϕ during the failed convergence was kept and fixed, and γ was then optimized along with the other parameters.

This procedure always allowed reaching convergence. It was preferred over an optimization with eight free parameters which would have been more difficult to realize: indeed, multiple optima were typically available with only seven free parameters as will be explained in the discussion part.

Carbon fiber microelectrodes modelling

Fitting of microelectrodes' cyclic voltammograms

The representation of cyclic voltammograms was done according to the method presented by

Compton, Laborda and Ward [START_REF] Compton | Understanding Voltammetry[END_REF]: the diffusion equation in cylindrical geometry (equation 11) was numerically solved using finite element methods with backward implicit method following Thomas' algorithm, as described in the authors' book.

(

) 11
Diffusion coefficient of reductant and oxidant were assumed to be of same value. An expanding grid was used for radial dimension. The first spatial boundary condition was equality of concentration with bulk concentration at maximum distance from the electrode:

(12)
The second spatial boundary condition was equality of concentration flux at the electrode's surface with consumption flux given by extended Butler-Volmer equation, as presented in equation ( 13):

The time boundary condition is given by equality of concentrations with bulk concentration at initial time.

(14)

Once the concentration profile is plotted along the spatial grid at different wanted times, the current profile can be deduced by the surface gradient according to equation ( 15):

The fitting of cyclic voltammetry experimental data was done as follows: experimental time array and measured potential array were taken as input variables, and corresponding response current was determined as a function of parameters to optimize. Then, typically the third cycle was used for the fitting: the calculated current points at corresponding input potentials were compared to experimental current points with least square minimization method. Optimized parameters were charge transfer kinetic constant k 0 , as well as geometric surface area A in order to include the measurement error of fibers' lengths as well as the variability of their radii. The whole optimization procedure was coded with MATLAB®.

Impedance model in cylindrical geometry for carbon fibers microelectrodes EIS

Carbon fibers microelectrodes can be reasonably represented as cylindrical electrodes and impedance at their surface can then simply be modelled with a Randles circuit. However, semi-infinite diffusion in cylindrical geometry leads to a Warburg impedance equation involving complicated Bessel functions [START_REF] Jacobsen | Diffusion impedance in planar, cylindrical and spherical symmetry[END_REF]. Landon-Lane and co-authors mentioned problems in convergence when using this expression and suggested the use of a simplified equation with errors in real and imaginary parts not exceeding 4% of their original values [START_REF] Landon-Lane | EIS at carbon fiber cylindrical microelectrodes[END_REF]. The total impedance, as written by aforementioned authors, is presented in equation ( 16)

and (17):

(16) (17)
Most symbols refer to the same physical parameters than the ones used for the macrohomogeneous model, and will therefore not be detailed here. The only exception is the CPE pseudocapacitance Y 0 .

In most cases, the assumption of a Nernstian system was verified (impedance minimum at estimated formal potential, and reversible-shaped cyclic voltammogram) and surface concentrations were obtained using Nernst equation with equations ( 8) and (9). However, in the case of non-activated fibres in (NH 4 ) 4 Fe(CN) 6 , that is for (NH 4 ) 4 Fe(CN) 6 with nonactivated fibers, the system did not appear Nernstian but quasi-reversible, based on the shape of the voltammograms and the impedance responses. Thus, kinetics became too slow for surface concentrations to be approximated with values given by Nernst equation, and mass transfer had to be considered to deduce their values.

In order to do so without introducing excessive complexity, steady-state and finite diffusion were assumed, semi-infinite diffusion assumption leading to no simple solution for diffusion equation in cylindrical geometry. After a development, presented in supplementary material, surface concentrations can be expressed as a function of mass transfer coefficient m O (assumed equal for oxidant and reductant):

(18) (19)
Parameters fitted were: R ohm , Y 0 , γ, k 0 and D. m O was also fitted in the case where surface concentrations could not be approximated to Nernst concentrations.

Results and discussion

Experimental results with porous "macro-felts"

4.1.1. Cyclic voltammetry on porous "macro-felts" Typical cyclic voltammograms on graphite felts are presented for both studied salts in

Figure 1, where the potential distance between oxidative and reductive peaks appear narrower with TEMPO-based molecule than with ferrocyanide. However, much dispersion was observed, especially with ferrocyanide, and distance between peaks could range from 20 to 100 mV with the same molecule (Figure S3). This dispersion is attributed to felts inhomogeneities. The height of the peaks, i.e. the peak currents, although similarly subject to some dispersion, was reproducibly different for the two molecules. Current peaks in TEMPObased molecules are about 1.3 times higher in intensity. On the other hand, a much higher non-faradaic current is observable on the CV plotted with activated felt, which corresponds to a higher amount of species adsorbed. This is well correlated by the increase in specific surface area measured by BET method: indeed, specific surfaces measured were about 0.5 m 2 g -1 for non-activated felt and 1.4 m 2 g -1 for activated one (about 500 cm 2 cm -3 and 1400 cm 2 cm -3 ).

Care should be given however to the meaning of these values, as these correspond to values referring to N 2 adsorption, and might represent a surface not readily available to bigger molecules used in this study.

Using Polarographica software, manual fits of the cyclic voltammograms of the two systems on non-activated felts are presented in Figure 2, estimated values for k 0 and total surface area A total are displayed below the graphs. The fit reproduces quite well the experimental data, especially around formal potential. Far from this value, important deviations are observed, which may arise from the non-faradaic current that is not considered in the model. Therefore, no fit was attempted for activated felts. Furthermore, the kinetic constants estimated here are actually a "minimum possible value" as the curves do not change shape anymore when inputting higher k 0 values. Total surface areas are different for the two systems and in agreement with observed differences in current peak values. The results of 45 cm 2 for ammonium hexacyanoferrate and 60 cm 2 for TEMPO derivative can be converted in specific surface areas using working electrode dimensions: a 1 cm diameter cylinder with length of 4.6 mm yields specific surface areas of 122 cm 2 cm -3 and 162 cm 2 cm -3 respectively. These values are of the same order of magnitude than the ones given by BET method, although sensibly lower, which could arise for example because of gases used with said method accessing extremely small pores on fibers surfaces. The total surfaces found here are also very close to the one found by other authors, although with felts from different providers [START_REF] Landon-Lane | Single fibre electrode measurements -A versatile strategy for assessing the non-uniform kinetics at carbon felt electrodes[END_REF], where 86.2 cm 2 was estimated for a 1 cm × 1 cm square piece of felt while in this study a surface between 45 cm 2 and 60 cm 2 corresponds to a 1 cm diameter cylinder. Nyquist diagrams show fast redox systems, conformingly to what could be observed with cyclic voltammograms: it is very hard to discern an impedance part corresponding to the "kinetic arc" that should typically be observed [START_REF] Lasia | Modeling of Impedance of Porous Electrodes[END_REF], [START_REF] Paasch | Theory of the electrochemical impedance of macrohomogeneous porous electrodes[END_REF], especially for activated felts.

Comparing the curves at other potentials around estimated formal potential allows to see more clearly that the area between the highest and intermediate frequencies ought to contain kinetic information (Figure S6).

Values of fitted parameters are available in Table 2 and all three replication's data is given in supplementary material. These values are presented under the form of the average of each fitted parameter calculated from each replication's fit, as well as the corresponding absolute standard deviation. Charge transfer kinetic constants k 0 obtained by the fit have values of the order of 10 -3 cm s -1 for both molecules on non-activated felt, with significantly higher values for ferrocyanide, and values of the order of 3.0-5.010 -2 cm s -1 for both molecules with activated felt. The relative deviation is high, between 33% and 50% of its value, except for TEMPO-based molecule on non-activated felt where it is very small. Values for non-activated felts are in good agreement with the order of magnitude found in the literature for similar materials [START_REF] Landon-Lane | EIS at carbon fiber cylindrical microelectrodes[END_REF], [START_REF] Landon-Lane | Single fibre electrode measurements -A versatile strategy for assessing the non-uniform kinetics at carbon felt electrodes[END_REF], [START_REF] Mccreery | Advanced carbon electrode materials for molecular electrochemistry[END_REF]. Some authors found extremely higher kinetic constant values in the order of 10 -1 cm.s -1 on the same felt without any mention of activation [START_REF] Luo | Unraveling pH dependent cycling stability of ferricyanide/ferrocyanide in redox flow batteries[END_REF], but this measurement was based on Nicholson analysis which assumes flat electrodes. Specific surface areas range between 150 and 350 cm 2 cm -3 , which is in good agreement with BET measurements on nonactivated felts and in even better agreement with estimation from manual fit of cyclic voltammograms. However, no increase in fitted surface area value is observed for activated felts, which does not agree with experimental BET results as well as cyclic voltammetry.

CPE exponent, when needed, would only slightly deviate from 1. Yet, this parameter was necessary for the convergence of kinetic constant to occur for activated felt. With such a high exponent value, typically above 0.97, capacitance value can be approximated by fitted pseudo-capacitance without the need of Brug relation. Its value is very high, ranging from 70 µF cm -2 to 100 µF cm -2 . Values of the order of 10 µF cm -2 or lower are expected for nonactivated graphite felt. Tortuosity coefficient f p was fitted, as its value was not known, and it also permitted to correct values of the felt porosity and the electrolyte resistivity. Its value was expected to be close to 1, but significantly higher values were needed to obtain a good fit at higher frequencies for non-activated felts. The parameters l p and ϕ were mainly considered adjustment parameters, much like f p , and no regard is given to their possible physical meaning. It is nevertheless possible to say that their values are very consistent: l p for the nonactivated felt present a value of about 10 µm whereas for activated one it shows a value of about 20 µm and the non-uniform diffusion coefficient ϕ is around 0.91 except for ferrocyanide on activated felt where it is about 0.82 with a higher deviation. We do not have any good explanation on why the slope of complex impedance at low frequencies is different for ferrocyanide/ferricyanide system on non-activated felt. This result was nevertheless systematically reproduced.

Therefore, Paasch and co-authors' model allows a very good fit of impedance data at the cost of many parameters to optimize. Consistent values for heterogeneous charge transfer kinetic constant and specific surface estimates are obtained but resulting double-layer capacitance is most probably overestimated. Other optima could have been selected, starting from different initial points in the optimization routine. Some optima were found to give a double-layer capacitance value for non-activated felt about ten times lower that the ones presented here, thus much closer to what should be expected. Optimized kinetic constant was then typically higher by one order of magnitude, and tortuosity coefficient closer to one, but the fit at lower frequencies was significantly worse. It came to mind that the ferrocyanide peculiarity at low frequency could actually mean a much slower kinetic constant, and it was indeed possible to find optima with kinetic constant of the order of 10 -4 cm s -1 for this system, but the specific surface area would then typically converge to a very high value that seemed absurd based on BET measurements and thus this result was rejected.

Experimental results with carbon fiber microelectrodes

Cyclic voltammetry on carbon fiber microelectrodes

Cyclic voltammograms on carbon fiber microelectrodes are presented for both studied salts in activated fiber exhibits a typically reversible curve, much like the curves obtained from TEMPO-based molecule, but non-activated fiber displays a slower system behavior. This result, corroborated with EIS measurements at different potentials, motivated the use of an adapted quasi-reversible model for the latter case.

In Figure 5 an example of experimental voltammogram with its fit is shown for each system on activated and non-activated fibers. The agreement between the fit and experimental data is overall very good, especially for activated fibers with fast kinetics, but some error can be seen typically for the reverse scan. Optimized parameters, k 0 and A, are presented in Table 3, as the average of fitted values for the three replications along with standard deviation. Kinetic constant values are indeed very low, of the order of 10 -5 cm s -1 , for ferrocyanide on nonactivated fibers. However, the activation permits a dramatic increase of charge transfer kinetic constant with this system, where kinetic constants reach several 10 -3 cm s -1 . Such a low kinetic constant with ferrocyanide/ferricyanide on carbon fiber microelectrodes was also observed by Landon-Lane, along with the dramatic increase after thermal activation [START_REF] Landon-Lane | Investigating Electron Transfer Kinetics in Redox Flow Batteries[END_REF]. The kinetic constants obtained for TEMPO-based molecule on non-activated fibers range from 310 -4 cm s -1 to almost 10 -3 cm s -1 , and after activation these values increase by an order of magnitude. The surface area was fitted to almost the same value as the geometric one in the case of TEMPO-based molecule, and for ammonium hexacyanoferrate in the case of activated fiber. However, its fitted value was slightly lower in the case of non-activated fibers for the latter system. This could not be the conjunction of an overestimated measurement of the fiber length and a fiber radius deviating from its average value, as these electrodes were the same as the ones for measurements with TEMPO-based system. As it was not the purpose of this study, this discrepancy was not further investigated.

Electrochemical impedance spectroscopy on carbon fiber microelectrodes

Nyquist and Bode phase diagrams for the two systems are presented in Figure 6. Here impedance was first normalized to the corresponding microelectrode geometric area for standard deviations to be meaningful. Otherwise, formalism used for the presentation of data in Nyquist diagrams is identical to the one previously described for graphite felts. The three replications' diagrams can be seen in Figure S7. The shape of all Nyquist diagrams is the one expected and reported in other works [START_REF] Landon-Lane | EIS at carbon fiber cylindrical microelectrodes[END_REF], [START_REF] Jacobsen | Diffusion impedance in planar, cylindrical and spherical symmetry[END_REF]: A half-circle containing charge transfer kinetics information is present and clearly identifiable at higher frequencies, while at lower frequencies complex impedance represents mass transfer and the imaginary part appears to tend towards a constant.

Optimized parameters' values, averaged over the three replications, are displayed in Table 4 and data for all replications is given in supplementary material. Optimized parameters with activated fibers all present significantly higher standard deviation values than non-activated felts. All k 0 obtained from EIS are in very good agreement with the ones obtained by fitting cyclic voltammograms: charge transfer kinetic constants k 0 is only slightly lower than 10 -3 cm.s -1 for TEMPO-based system with non-activated fibers and of the order of several 10 -3 cm.s -1 with activated ones, while the diffusion coefficient D is in both cases about 410 -6 cm 2 .s -1 which is basically the value that was determined with other methods (results not shown in this work). For ferrocyanide/ferricyanide system, k 0 is of the order of 10 -5 cm.s -1 with non-activated fibers whereas with activated fibers it has about the same value and standard deviation as the TEMPO-based molecule, which makes for a dramatic increase of a factor 500 due to activation, as was observed when fitting cyclic voltammograms. This validates the specific formalism that was used to fit the spectra for ammonium ferrocyanide.

The diffusion coefficient is optimized at a value close to its measured one for activated fibers, but is far from the correct value of 4.510 -6 cm 2 .s -1 in the case of non-activated fibers. It was attempted to fix this value, but resulting optima would yield poor fitting of experimental data.

On the other hand, the mass transfer coefficient fitted with quasi-reversible model in the case of (NH 4 ) 4 Fe(CN) 6 with non-activated fiber is 1.010 -3 cm.s -1 , and this value is not absurd when compared to typical values of this parameter [START_REF] Bard | Electrochemical Methods[END_REF], [START_REF] Gao | Rotating Microdisk Voltammetry[END_REF]. Furthermore, using Szabo's approximation for semi-infinite diffusion in cylindrical geometry [START_REF] Szabo | Chronoamperometric current at hemicylinder and band microelectrodes: Theory and experiment[END_REF] in the limiting case of quasi-steady-state as reported by Bard and Faulkner [31], and as presented in equation ( 20), it is possible to estimate a mass transfer coefficient of about 2.910 -3 cm.s -1 taking a value of five minutes for t (buffer time applied before each potential step in EIS). Taking a time of one hour only changes the value to 2.310 -3 cm.s -1 and these values agree very well with the one optimized here.

The fitted values of CPE element's parameters are similar for the two redox systems studied here: the exponent γ is very close to 1 for non-activated fiber but decreases to a value of about 0.8 for activated fibers. The deviation of this parameter from unity is usually attributed to a more inhomogeneous electrode surface [START_REF] Brug | The analysis of electrode impedances complicated by the presence of a constant phase element[END_REF], [35], which is expected for a felt having undergone for instance a thermal treatment. Non-activated fibers have a pseudo-capacitance of several µF.cm -2 , which is very small but not unexpected for this material [START_REF] Kovach | Electrochemistry at partially blocked carbon-fiber microcylinder electrodes[END_REF], [START_REF] Miller | Kinetic Study of Electrochemical Treatment of Carbon Fiber Microelectrodes Leading to In Situ Enhancement of Vanadium Flow Battery Efficiency[END_REF], [START_REF] Landon-Lane | EIS at carbon fiber cylindrical microelectrodes[END_REF].

Activated ones show dramatically higher values for this parameter, with about 25 µF.cm -2 for TEMPO-based system and going above 100 µF.cm -2 for ammonium hexacyanoferrate. One may recombine Brug relation [START_REF] Brug | The analysis of electrode impedances complicated by the presence of a constant phase element[END_REF] to access actual capacitance values from those estimated:

Capacitances on activated felts obtained with this method are between 6 and 25 µF.cm -2 for the TEMPO-based compound and between 40 and 130 µF.cm -2 for ferrocyanide. This increase in double-layer capacitance represents a higher number of charges adsorbed at the electrode's surface, which could arise from an increase of this surface. This is also observed and described in Landon-Lane's work [START_REF] Landon-Lane | Investigating Electron Transfer Kinetics in Redox Flow Batteries[END_REF].

Comparison between macro-felts and micro-fibers

The charge transfer kinetics of two redox systems on activated and non-activated carbon felts are studied at two different scales: macro-felt and micro-fiber. For each scale, CV and EIS are carried out and analyzed by corresponding models. At macro-felts, CV analysis is based on Polarographica software and EIS data are fitted with a modified macrohomogeneous model.

At micro-fibers, CV is fitted by considering non-Nernstian charge transfer kinetics and 1D semi-infinite diffusion in cylindrical coordinates, and EIS is fitted with a modified Randles circuit adapted to cylindrical geometry.

For macro-felts, the oxidation and reduction peak separation is an indication of charge transfer kinetics on the electrode. From the results in Figure 1, it is seen that the differences between activated and non-activated felts in both TEMPO-based compound and ferrocyanide electrolytes are insignificant. This is confirmed by Polarographica simulations where voltammograms almost do not change when increasing k 0 to values higher than 10 -2 cm.s -1 .

EIS is more sensitive to differences in charge transfer kinetics. Qualitatively, differences can be observed from the Nyquist plot in high-frequency range (Figure 3). By fitting with the macrohomogenous model, k 0 value is one order of magnitude higher for activated felts than for non-activated ones in both TEMPO-based and ferrocyanide electrolytes and this is consistent with what can be found for other systems in literature [START_REF] Mazúr | Performance evaluation of thermally treated graphite felt electrodes for vanadium redox flow battery and their four-point single cell characterization[END_REF].

For micro-fibers, the difference in kinetics yields a clear change in the shape of CV, especially for ferrocyanide electrolyte. The k 0 values derived from fitting voltammograms show an increase of 300 and 10 times by activation in ferrocyanide and TEMPO electrolytes, respectively. Furthermore, they are in good agreement with the values analyzed from EIS measurements (Table 4). This suggests that both CV and EIS may provide reliable kinetic information of the micro-fibers. We consider that the micro-fibers have two main advantages over the macro-felts for the quantitative analysis: well-defined geometry that gets rid of the approximations in felt structure and simplifies the mathematical treatment in 1D, and much lower capacitive contribution in CV due to much lower active surface area. However, it must be pointed out that k 0 values analyzed from micro-fibers are systematically lower than that analyzed from macro-felts where the differences may reach ca. two orders of magnitude in some cases. The reason for these discrepancies is unclear but one possible explanation could be the analytical model for macro-felts, where the geometrical approximations might not always correspond well to the experimental systems. From the practical point of view, direct measurements on macro-felts are more convenient and straightforward. On the other hand, measurements on extracted micro-fibers may reveal more sensitively the kinetic difference for electrode materials undergoing different treatments. Moreover, they may also provide information on the heterogeneity of the felt by extracting fibers randomly from different areas.

Conclusion

Charge transfer kinetic constant was measured using CV fitting or simulation and EIS on both macro-felts and micro-fibers, activated or not, for two aqueous organic electrolytes: ferrocyanide and TEMPO-based molecule. For macro-felts, EIS is more sensitive to k 0 determination than CV analysis in the range of k 0 considered, and an increase by one order of magnitude is measured for activated material. For micro-fibers, differences can already be observed on CV curves, and analysis using CV fitting and EIS give k 0 values with very good consistency. However, k 0 measured on micro-fibers are typically lower than the ones obtained on macro-felts. The factor by which these values are decreased is usually about one order of magnitude, but reaches about two orders of magnitude in the case of ferrocyanide on nonactivated felt. Since results given by the two electrochemical methods are consistent with each other for the felt and for the micro-fiber, the origin of this discrepancy is not well understood.

Since physical analysis on microelectrodes is much straightforward and subject to less approximations, it is reasonable to attribute said discrepancy to possible misrepresentation of the complex porous structure by the models used. Because of the high sensitivity of methods used on micro-fibers to detect k 0 variation with material treatment, this support makes for an excellent tool to detect the effect of various treatments on charge transfer kinetics for carbonbased felt materials.

[s -1 ] Imaginary number j [-] CPE exponent γ [-] Non-uniform diffusion exponent ϕ [-] Ohmic resistance R ohm [Ω] Charge transfer resistance R ct [Ω] Carbon fiber microelectrode radius r 0 [cm] Pseudo-capacitance (for microelectrode model) Y 0 [F cm -2 s (γ -1) ] Mass transfer coefficient m O [cm s -1 ]
Figures:

(a Tables: Y 0 [F cm -2 s γ-1 ] 2.6810 -6 5.1010 -7 2.5610 -5 1.5810 -5 3.6410 -6 2.0510 -6 1.0810 -4 5.6310 k 0 [cm s -1 ] 9.0010 -4 1.7510 -4 5.8110 -3 3.6310 -3 1.7610 -5 6.4410 -6 5.9510 -3 3.8410 -3 D [cm 2 s -1 ] 4.4710 -6 1.4410 -6 3.8810 -6 5.5910 -7 1.2610 -5 5.5810 -6 4.3810 -6 9.1710 -7 m O [cm s -1 ] ----1.0010 -3 1.6310 -4 --

  C O s and C R s are the surface concentrations of respectively the oxidant and the reductant of considered redox couple; studied salts were assumed to have a Nernstian behavior on graphite felts, and therefore these values were obtained through Nernst equation as follows:

4. 1 . 2 .

 12 Electrochemical impedance spectroscopy on porous "macro-felts" Nyquist and Bode phase diagrams at formal potential are displayed in Figure 3 for TEMPObased and ferrocyanide molecules. Each experimental diagram corresponds to the replication giving intermediate values, out of the three performed, along with its fit using macrohomogeneous model. The three replications' diagrams can be seen in Figure S4. A strong dispersion on dc resistance value was observed, which resulted in Nyquist diagrams having very close shapes despite being sometimes far from one another on the real numbers' axis, which is proven in Figure S5 by translating the curves on the horizontal direction.

Figure 4 .

 4 Figure 4. While the curves for TEMPO-based molecule are close in shape for non-activated

1 :

 1 ) TEMPO-based molecule (b) Ammonium hexacyanoferrate Figure Typical cyclic voltammograms on SGL graphite felts with (a) TEMPO-based molecule and (b) ammonium hexacyanoferrate. Scan rate: 10 mV.s -1 (a) TEMPO-based k 0 = 1.0 10 -2 cm.s -1 A total = 60 cm 2 (b) (NH 4 ) 4 Fe(CN) 6 k 0 = 1.0 10 -2 cm.s -1 A total = 45 cm 2

Figure 2 : 3 : 4 : 5 : 6 :

 23456 Figure 2: Simulation of cyclic voltammetry with "Statistically Weighted External Cylindrical Finite Diffusion" model using Polarographica software (in red) compared with compensated experimental data (in black) on nonactivated SGL graphite felt with TEMPO-based system (left) and ammonium hexacyanoferrate (right).Scan rate: 10 mV.s -1 . Concentration: 10 -4 M in supporting electrolyte.

  10.1016/j.jpowsour.2018.01.079. [37] R. C. Weast and D. R. Lide, CRC Handbook of Chemistry, and Physics, 70th ed. 1989.

	Appendix: List of symbols
	Parameter	Symbol Unit
	Number of electrons exchanged	n [-]
	Faraday constant	F [C mol -1 ]
	Ideal gas constant	R g [J mol -1 K -1 ]
	Temperature	T [K]
	Charge transfer coefficient	α [-]
	Oxidant concentration near the electrode surface Reductant concentration near the electrode surface Oxidant concentration in bulk electrolyte Reductant concentration in bulk electrolyte Formal potential of the redox system	s [mol cm -3 ] C O s [mol cm -3 ] C R * [mol cm -3 ] C O * [mol cm -3 ] C R E 0 ' [V]
	dc potential applied before EIS measurement	Ē [V]
	Carbon felt porosity	ε [-]
	Carbon felt resistivity Bulk electrolyte resistivity Macrohomogeneous model parameter 1	ρ 1 [Ω cm] * [Ω cm] ρ 2 k [s -1 ]
	Macrohomogeneous model parameter 2	K [cm 2 s -1 ]
	Macrohomogeneous model parameter 3	 [-]
	Hindrance function	y [-]
	Double-layer capacitance	C [F cm -2 ]
	Specific surface area Heterogeneous charge transfer kinetic constant	S c [cm 2 cm -3 ] k 0 [cm s -1 ]
	Diffusion coefficient	D [cm 2 s -1 ]
	Diffusion layer thickness	l p [cm]
	Porous electrode thickness	d [cm]
	Electrode geometric area (cross section for carbon felt, cylinder area for microelectrode)	A [cm 2 ]
	Time	t [s]
	Angular frequency	ω

Table 1 :

 1 Physical parameters fixed for the macrohomogeneous model Gross estimation based on literature conductivity values of support electrolyte [37]f measured from Linear Scan Voltammetry on a glassy carbon rotating electrode (not shown in this work)

	Parameter	Symbol Unit	TEMPO-based	(NH 4 ) 4 Fe(CN) 6
	Charge transfer coefficient	α [-]	0.5	0.5
	Carbon felt Porosity	ε [-]	0.95 a	0.95 a
	Carbon felt resistivity Bulk electrolyte resistivity Diffusion coefficient	ρ 1 [Ω cm] * [Ω cm] ρ 2 D [cm 2 s -1 ]	0.5 b 12.5 c 3.8 10 -6 f	0.5 b 10 d 4.5 10 -6 f
	a,b taken from providers' estimations
	c,d

Table 2 :

 2 Averaged optimized parameters' values (calculated from the three replications) and their standard deviation values for EIS experiments carried out on SGL graphite felt with TEMPO-based molecule (left) at estimated formal potential of 0.66 V vs Ag/AgCl and with ammonium hexacyanoferrate (right) at estimated formal potential of 0.25 V vs Ag/AgCl Average Std Dev Average Std Dev Average Std DevC [F cm -2 s γ-1 ] 1.0210 -4 3.5110 -6 9.5610 -5 1.3410 -5 5.2710 -5 7.3410 -6 9.0410 -5 2.2110-6

			TEMPO-based			(NH) 4 Fe(CN) 6
		Non-activated	Activated	Non-activated	Activated
	Average Std Dev R ohm [Ω] 19.1 17.7	31.1	2.3	64.0	48.2	51.2	41.1
	γ [-]	0.97	0.002	0.99	0.01	1.00 *	0.00	1.00	0.00
	k 0 [cm s -1 ]	1.5710 -3 7.2610 -5 5.7710 -2 3.9710 -2 2.2010 -3 3.7010 -4 3.2710 -2 2.2610 -2
	l p [cm]	0.0011	0.00003	0.0020	0.0001	0.0010	0.00009	0.0019	0.00004
	S c [cm 2 cm -3 ]	331.8	4.1	286.4	23.4	294.6	9.6	278.1	35.7
	f p [-]	3.01	0.24	1.21	0.11	6.44	0.87	1.31	0.21
	ϕ [-]	0.91 *	0.00	0.91 *	0.01	0.81	0.02	0.91 *	0.01

* value fixed during optimization

Table 3 :

 3 Average of optimized parameters' values (calculated from the three replications) and their standard deviation values for microelectrodes' cyclic voltammograms with TEMPO-based compound (left) and ammonium hexacyanoferrate (right). The resulting ratio between optimized surface area and measuredgeometric area is also given.

			TEMPO-based			(NH) 4 Fe(CN) 6
		Non-activated	Activated	Non-activated	Activated
		Average Std Dev Average Std Dev	Average Std Dev Average Std Dev
	k 0 [cm s -1 ] 6.8710 -4 3.2210 -4 6.5810 -3 4.1610 -3 1.8910 -5 9.3210 -6 4.7610 -3 2.4610 -3
	A [cm 2 ]	7.1010 -4 1.4110 -4 8.0810 -4 2.0310 -4 6.4210 -4 8.8210 -5 7.7910 -4 1.3210 -4
	A fit /A geo [-]	0.96	0.14	1.07	0.05	0.85	0.19	1.05	0.12

Table 4 :

 4 Averaged optimized parameters' values (calculated from the three replications) and their standard deviation values for EIS experiments carried out on carbon fiber microelectrode drawn from SGL graphite felt with TEMPO-based molecule (left) at estimated formal potential of 0.66 V vs Ag/AgCl, and for ammonium hexacyanoferrate (right) at 0.60 V vs Ag/AgCl for non-activated fibers and at estimated formal potential of 0.25 V vs Ag/AgCl for activated fibers. Average Std Dev Average Std Dev Average Std Dev

	TEMPO-based			(NH) 4 Fe(CN) 6
	Non-activated	Activated	Non-activated	Activated
	Average Std Dev R ohm [Ω] 534.4 172.2	2605.0	3736.3	569.0	187.4	2090.4	1555.5

Acknowledgements

This work was carried out under Higreew project, which has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement no. 875613. The authors are grateful to Dr. Tim Tichter and M.Sc. Jonathan Schneider from Physikalische und Theoretische Chemie, Institut für Chemie und Biochemie, Berlin, Germany, for insightful discussion as well as for providing us with their excellent software: Polarographica. The authors wish to thank M.Sc. Oihane Zugazua and Drs. Maddalen Agirre, Ana Catarina Lopes and Eduardo Sánchez-Díez, from CICEnergigune, for providing the electrolytes studied in this paper.


Josh A Taylor
email: josh.taylor@utoronto.ca

Alain Rapaport
email: alain.rapaport@inrae.fr

Denis Dochain
email: denis.dochain@uclouvain.be.

Convex optimization of bioprocesses

Keywords: Bioprocess, compartmental system, convex relaxation, second-order cone programming, wastewater treatment

We optimize a general model of bioprocesses, which is nonconvex due to the microbial growth in the biochemical reactors. We formulate a convex relaxation and give conditions guaranteeing its exactness in both the transient and steady state cases. When the growth kinetics are modeled by the Contois or, under constant biomass, Monod or Powell functions, the relaxation is a second-order cone program, which can be solved efficiently at large scales. We implement the model on a numerical example based on a wastewater treatment system.

I. INTRODUCTION

We optimize a model of dynamical bioprocesses consisting of a set of biochemical reactors interconnected by diffusion and mass flow. The objectives include minimizing substrate outflow, maximizing biogas production, and tracking setpoints. Within each reactor, several microbial reactions convert any number of biotic or abiotic reactants into biomass and/or products. This setup describes a variety of physical systems such as wastewater treatment networks, the production of various chemicals, and compartmental approximations of bioprocesses in continuous media.

The optimization is nonconvex due to the nonlinear microbial growth in the reactors. We formulate a convex relaxation and derive analytical conditions guaranteeing its exactness, which is to say the relaxed solution also solves the original problem. We focus on the case when the microbial growth can be represented as a second-order cone (SOC) constraint, as recently shown in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF] for the Monod [START_REF] Monod | The growth of bacterial cultures[END_REF] and Contois [START_REF] Contois | Kinetics of bacterial growth: relationship between population density and specific growth rate of continuous cultures[END_REF] growth rates. We further show that the Powell function [START_REF] Powell | The growth rate of microorganisms as a function of substrate concentration[END_REF] also admits an SOC representation when the biomass is constant. The resulting optimization is a second-order cone program (SOCP), which can be solved at large scales with commercial software.

The most closely related topics to ours are chemical process optimization, control and optimization of wastewater systems, and control of bioprocesses. Most existing approaches to process optimization [START_REF] Biegler | Nonlinear programming: concepts, algorithms, and applications to chemical processes[END_REF] and wastewater [START_REF] Ocampo-Martinez | Model predictive control of wastewater systems[END_REF] do not explicitly model the microbial growth, and often use either linear programming or general nonlinear solvers. There have been many applications of nonlinear control [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] and optimization [START_REF] Srinivasan | Dynamic optimization of batch processes: I. characterization of the nominal solution[END_REF] to bioprocesses, but not convex relaxations or SOCP.

Our main results are generalizations of those in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], which focused on the gradostat with a single reaction [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]. Here we allow for any number of substrates and biomasses, general convex objectives, and multiple biochemical reactions. Our original theoretical contributions are as follows.

• In Section III, we formulate convex relaxations for optimizing the trajectory and steady state solution of a general bioprocess.

• In Section III-C, we identify an SOC representation of the relaxed growth constraint when the growth rate is the Powell function [START_REF] Powell | The growth rate of microorganisms as a function of substrate concentration[END_REF], [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] and the biomass is constant. • In Section IV, we give conditions under which the relaxations are guaranteed to be exact in both the transient and steady state cases.

To streamline exposition, the only external inputs to the model are the influent concentrations, e.g., biochemical oxygen demand and ammonia. Our main exactness results straightforwardly apply when the flow rates are also variable, but the resulting bilinearities make the problem nonconvex. This can be handled using techniques like disjunctive programming, as in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], or further convex relaxation, which we discuss at the end of Section II-E.

We apply our results in two examples. In Section V-A, we show that our exactness conditions simplify to those in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF] when specialized to the gradostat. In Section V-B, we optimize the allocation of sewage to three wastewater treatment plants over two weeks. The relaxation is exact, and takes roughly twenty minutes to solve using SOCP [START_REF] Lobo | Applications of second-order cone programming[END_REF].

II. SETUP

A. Network modeling

The system consists of s well-mixed tanks interconnected by mass flow and diffusion. We denote the set of tanks S. V ∈ R s×s is a diagonal matrix in which V ii is the volume of tank i. Tank i has water inlet flow rate Q in i and outlet flow rate Q out i . We let Q ij denote the flow from tank i to tank j. Let d ij denote the diffusion rate between tanks i and j, where

d ij = d ji . Let C = diag Q in i , M ij = Q ji , i = j -Q out i -k∈S Q ik , i = j , L ij = d ij , i = j -k∈S d ik , i = j ,
and N = M + L. M and L are respectively compartmental and Laplacian matrices. M is invertible if the network is outflow connected, which is to say that there is a directed path from every tank to some tank with outflow [START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF]. Because L is negative semidefinite, N is also invertible if M is outflow connected, and potentially even if M is not outflow connected.

B. Microbial growth

We model the microbial growth in the tanks using the notation of Section 1.5 of [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]. There are m substrates and biomasses in each perfectly mixed tank. ξ i ∈ R m + is the process state vector of tank i ∈ S, which contains the concentrations of the substrates and biomasses, and ξ in i ∈ R m + is the corresponding influent concentration vector. This model is minimal in that ξ includes intermediary products, e.g., substrates produced by one reaction and consumed by another, but not final products such as the CH 4 ultimately produced by anaerobic digestion.

There are r different types of reactions that convert substrates to other substrates and biomasses. φ i (ξ i ) ∈ R r + is a vector of the reaction kinetics in tank i. We are interested in the case where the elements of φ i (ξ i ) are concave functions, and in particular representable as SOC constraints. We show how to do this for Monod and Contois kinetics later in Examples III-A and III-B.

Let κ i ∈ R m×r be the stoichiometric matrix relating the reaction vector, φ i (ξ i ), to the evolution of the process state in tank i. The dynamics in tank i ∈ S are

V ii ξi = V ii κ i φ i (ξ i ) -Q out i ξ i - j∈S (Q ij + d ij )ξ i + Q in i ξ in i + j∈S (Q ji + d ij )ξ j .
The following example illustrates ξ i and κ i .

Example 1 (Two-step anaerobic digestion): There are two substrates, S a i and S b i , and two biomasses, X a i and X b i . We let

ξ i = S a i , S b i , X a i , X b i . S a i is converted to both X a i and S b i at the rate µ a (S a i ) X a i . S b i is converted to X b i at the rate µ b S b i X b i . Therefore, φ i (ξ i ) = µ a (S a i ) X a i µ b S b i X b i and κ i =     -1 0 1 -1 1 0 0 1     .
We now write the dynamics in vector form. We suppress subscripts to represent stacked vectors, i.e., ξ = [ξ 1 , ..., ξs] and φ(ξ) = [φ 1 (ξ 1 ), ..., φs(ξs)] . Let A ⊗ B denote the Kronecker product of A and B, Iα ∈ R α×α the identity matrix, and  = A ⊗ Im. Let K be a block diagonal matrix with κ 1 , ..., κs on its main diagonal. If κ i = κ for all i ∈ S, then K = Is ⊗ κ. The dynamics of the full system are given in vector form by

V ξ = V Kφ(ξ) + N ξ + Ĉξ in . (1)
We allow the dynamics to be non-autonomous, in which case N , Ĉ, and ξ in can be time-varying.

C. Discretization in time

To make (1) compatible with finite-dimensional optimization, we replace the derivatives with a numerical approximation, which we denote Dn. For example, in the case of the implicit Euler method with time step ∆, Dn[ξ(•)] = (ξ(n) -ξ(n -1))/∆. Dn could also be a more sophisticated approximation such as a Runge-Kutta scheme [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF]. The time periods are indexed n ∈ N = {1, ..., τ }. We have

V Dn[ξ(•)] = V Kφ(ξ(n)) + N (n)ξ(n) + Ĉ(n)ξ in (n) (2)
for n ∈ N . The initial condition is ξ(0) = ξ 0 .

D. Objectives

We consider objectives of the form

F(ξ, T ) = n∈N F ξ (ξ(n)) + F φ (T (n)), (3)
where F ξ and F φ are convex and T (n) = φ(ξ(n)). The following are examples.

• Minimizing the outflow of substrates,

F ξ (ξ(n)) = i∈S Q out i η i ξ i (n),
where η i is a vector that selects the entries of ξ i (n) corresponding to pollutants. • The production of biogas in a tank is proportional to the kinetics that convert substrates to biomass. Let σ i ∈ R m + be a vector that is only nonzero for entries of T i (n) corresponding to biogas production. Let M ⊆ S be the subset of tanks that can capture biogas from anaerobic digestion. We maximize biogas through the objective

F φ (T (n)) = - i∈M V ii σ i T i (n). • Setpoint tracking, F ξ (ξ(n)) = ξ(n) -ξ A ξ(n) -ξ ,
where A 0 and ξ is a desired operating point.

E. Problem statement

We aim to solve the following optimization problem.

P min ξ,ξ in ,T F(ξ, T ) (4a) such that T (n) = φ(ξ(n)), n ∈ N (4b) V Dn[ξ(•)] = V KT (n) + N (n)ξ(n) + Ĉ(n)ξ in (n), n ∈ N (4c) ξ, ξ in , T ∈ Ω. (4d)
P models the optimization of a broad range of bioprocesses such as wastewater treatment. We refer the reader to [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] for broad coverage of this topic.

A solution of P is a trajectory ξ(n), ξ in (n), T (n) , n ∈ N . The flows and diffusions between tanks, encoded by the matrices N (n) and Ĉ(n), are not decision variables. For this reason constraint (4c) is linear. The set Ω in (4d) consists of linear constraints such as the initial condition, total input matter, and maximum substrate concentrations; several other examples are given in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF]. Note that Ω can constrain T (•) so as to allow constraints on the growth without adding nonlinearities. The only nonconvexity is therefore (4b), the growth constraint; this is the focus of the next two sections.

We note that in many applications, the flow rates are important decision variables. They are parameters here because our focus is on incorporating the growth kinetics in a convex fashion. In the case that the flow rates are variable, (4c) becomes bilinear and hence nonconvex. There are several ways to handle the bilinearity, including

• McCormick [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I-Convex underestimating problems[END_REF] and lift-and-project [START_REF] Sherali | A new reformulation-linearization technique for bilinear programming problems[END_REF] relaxations;

• disjunctive programming reformulations if the flow variables are binary, as in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF]; • and finding a local minimum via nonlinear programming, e.g., exploiting the biconvex structure with the Alternating Direction Method of Multipliers [START_REF] Boyd | Distributed Optimization and Statistical Learning Via the Alternating Direction Method of Multipliers[END_REF]. All three of the above techniques are viable because, as described in the next section, we have a tractable way to represent the growth constraint, (4b).

III. CONVEX RELAXATION

P is nonconvex because constraint (4b) is a nonlinear equality. One way around this difficulty is to instead solve a convex relaxation of P, as in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF]. If all elements of the vector φ(•) are concave functions, we obtain a convex relaxation by replacing (4b) with the inequality

T (n) ≤ φ(ξ(n)), n ∈ N . (5)
We refer to the resulting optimization as P R . As mentioned earlier, we are interested in the case where ( 5) is concave and, ideally, representable as linear or SOC constraints. Several such examples are given below. We refer the reader to Appendix 1 of [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] for a comprehensive list of growth rates.

A. Contois growth

Suppose that there is a substrate of concentration S, a biomass of concentration X, and the growth rate is Contois [START_REF] Contois | Kinetics of bacterial growth: relationship between population density and specific growth rate of continuous cultures[END_REF]. Then constraint (5) takes the form

T (n) ≤ µS(n)X(n) k C X(n) + S(n) ,
where µ and k C are constant parameters. As shown in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], this is concave and can be written as the SOC constraint

  µS(n) k C T (n) µk C X(n)   ≤ µk C X(n) + µS(n) -k C T (n). (6)

B. Monod growth with constant biomass

Consider Section III-A, but now suppose that the growth rate is Monod [START_REF] Monod | The growth of bacterial cultures[END_REF]. Then constraint (5) takes the form

T (n) ≤ µS(n)X(n) k M + S(n) .
This constraint is quasiconcave. It becomes concave if we assume that the biomass in each time period is not an optimization variable, but an exogenous parameter, i.e., X(n) = X(n) for n ∈ N . This approximation is often valid because the biomass concentration is typically larger and varies more slowly than the substrate concentrations, and is therefore relatively insensitive to the substrates. In this case, as shown in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], it can be written as the SOC constraint

  µS(n) X(n) k M T (n) µk M X(n)   ≤ µk M X(n) + µS(n) X(n) -k M T (n). (7)

C. Powell growth with constant biomass

In the case of the Powell growth rate [START_REF] Powell | The growth rate of microorganisms as a function of substrate concentration[END_REF], constraint (5) takes the form

T (n) ≤ µX(n) 2K B K B + S(n) -K B + S(n) 2 -4K B S(n) , (8)
where 0 < K B ≤ 1/4, with the upper limit ensuring the square root is real. The following is an original result.

Theorem 1: If X(n) = X(n) for n ∈ N , (8) is equivalent to the SOC constraint S(n) -2K B K B -4K 2 B ≤ S(n) + K B - 2K B µ X(n) T (n). ( 9
)
This can be shown by simplifying [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF].

D. Blackman growth with constant biomass

In the case of the Blackman growth rate [START_REF] Blackman | Optima and limiting factors[END_REF], constraint (5) takes the form

T (n) ≤ µ K B S(n)X(n) if S(n) ≤ K B µX(n) if S(n) > K B ,
where K B > 0. This is a pair of linear constraints if, as above, X(n) = X(n) for n ∈ N .

E. Interactive and non-interactive growth

Suppose the growth of the biomass depends on two rates, µ a (S(n), X(n)) and µ b (S(n), X(n)), and the individual kinetics constraints,

T a (n) ≤ µ a (S(n), X(n))X(n) and T a (n) ≤ µ b (S(n), X(n))X(n), both have SOC representations.
The dependency often takes one of two forms: noninteractive, min{µ a (S(n), X(n)), µ b (S(n), X(n))}, which in ecological modeling is known as Liebig's Law, and interactive, µ a (S(n), X(n))µ b (S(n), X(n)), which is common in models of bioprocesses.

The non-interactive case is enforced by the two individual kinetics constraints along with T (n) ≤ T a (n) and T (n) ≤ T b (n). The interactive case is in general nonconvex, and does not have an SOC representation. However, the geometric mean of the growth rates,

µ a (S(n), X(n))µ b (S(n), X(n))
, does lead to an SOC representation. It is enforced by the two individual kinetics constraints along with T (n) 2 ≤ T a (n)T b (n). The latter is hyperbolic, a type of SOC constraint [START_REF] Lobo | Applications of second-order cone programming[END_REF]. We note that if both rates are Monod functions, this is somewhat similar to the growth rate given in [START_REF] Jost | Interactions of Tetrahymena pyriformis, Escherichia coli, Azotobacter vinelandii, and glucose in a minimal medium[END_REF], which does not have a square root.

IV. EXACTNESS

P R is exact if its optimal solution also solves P; this is the ideal outcome. We can verify exactness by checking whether ( 5) is satisfied with equality at the optimal solution. When ( 5) is not satisfied with equality, P R might still provide a close approximation of P, but this is hard to guarantee. While the exactness of a solution can always be tested, it is also useful to have analytical conditions, even if narrow, under which the exactness of P R is guaranteed.

For the rest of this section we let

Dn[ξ(•)] = (ξ(n) -ξ(n -1))/∆,
which corresponds to the implicit Euler step. We assume that the only constraint specified by Ω in (4d) is an initial condition, ξ(0) = ξ 0 , and that ξ in (•) is fixed; note that as long as exactness holds for all feasible values of ξ in (•), it holds when ξ in (•) is a variable. We also assume that strong duality holds for P R ; this is a mild assumption that, e.g., holds as long as there is a feasible solution in which ξ(n) > 0 for all n ∈ N . Let J (ξ(n)) ∈ R rs×ms denote the Jacobian matrix of φ(•) at ξ(n). For convenience, we define the following quantities:

Γ(n) = 1 ∆ V /∆ -N (n) -J (ξ(n)) K V -1 V = Ims -∆ V -1 N (n) + J (ξ(n)) K V -1 Φ(n) = -∇F φ (T (n)) -∆K V τ k=n k l=n Γ(l) V -1 × ∇F ξ (ξ(k)) + J (ξ(k)) ∇F φ (T (k)) .
Observe that if ∆ is small enough, Γ(n) is positive definite and close to the identify matrix. Theorem 2: P R is exact if at an optimal solution, Φ(n) > 0 for all n ∈ N .

Proof: Let λ(n) ∈ R ms and ρ(n) ∈ R rs be the respective dual multipliers of constraints (4c) and [START_REF] Biegler | Nonlinear programming: concepts, algorithms, and applications to chemical processes[END_REF] for n ∈ N . The Lagrangian of P R is

L = F(ξ, T ) + n∈N ρ(n) (T (n) -φ(ξ(n))) + λ(n) V (ξ(n -1) -ξ(n))/∆ + V KT (n) + N (n)ξ(n) + Ĉ(n)ξ in (n) .
Differentiating the Lagrangian by T (n) and ξ(n) and setting it to zero gives

-ρ(n) = ∇F φ (T (n)) + K V λ(n), n ∈ N , (10a)
J (ξ(n)) ρ(n) = ∇F ξ (ξ(n)) -V /∆ -N (n) λ(n) + V λ(n + 1)/∆, n ∈ N \ τ (10b) J (ξ(τ )) ρ(τ ) = ∇F ξ (ξ(τ )) -V /∆ -N (n) λ(τ ). ( 10c
)
We now solve for ρ(n). Premultiplying (10a) by J (ξ(n)) and summing with (10b) and (10c) gives

λ(τ ) = ∆Γ(τ ) V -1 ∇F ξ (ξ(τ )) + J (ξ(τ )) ∇F φ (T (τ )) ,
and, for n ∈ N \ τ ,

λ(n) = ∆Γ(n) V -1 ∇F ξ (ξ(n)) + J (ξ(n)) ∇F φ (T (n)) + Γ(n)λ(n + 1).
Expanding the recursion yields

λ(n) = ∆ τ k=n k l=n Γ(l) V -1 × ∇F ξ (ξ(k)) + J (ξ(k)) ∇F φ (T (k))
.

We now substitute this into (10a) to obtain

ρ(n) = -∇F φ (T (n)) -∆K V τ k=n k l=n Γ(l) V -1 × ∇F ξ (ξ(k)) + J (ξ(k)) ∇F φ (T (k)) = Φ(n).
From here we can see that the conditions of the theorem guarantee that ρ(n) > 0 for all n ∈ N . Theorem 2 is of limited immediate use because we must know the optimal solution of P R to test if Φ(n) > 0. It can however be used to derive sufficient conditions for exactness that are easy to test. We now derive two such conditions that do not require knowledge of the optimal solution.

For the rest of this section, assume that F ξ and F φ are linear with gradients f ξ ∈ R ms and f φ ∈ R rs . In this case, we can write

Φ(n) = -∆K V τ k=n k l=n Γ(l) V -1 f ξ - Ims + ∆K V τ k=n k l=n Γ(l) V -1 J (ξ(k)) f φ .
We also assume that each element of the vector of reaction kinetics has bounded slope, so that all entries of J (ξ(n)), n ∈ N , are bounded.

Let

Ψ(n) = V -1 N (n) + V -1 J (ξ(n)) K V .
The first term of Ψ(n) is negative semidefinite because N (n) is compartmental [START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF]. The latter term is bounded by assumption, and as we will see in the examples, usually negative semidefinite-we assume that this is the case. We therefore assume that the eigenvalues of Ψ(n) are in the range [-ψ, 0], where ψ > 0 is an upper bound on the magnitude. We can use the push-through identity to write

Γ(n) = Ims + ∆Ψ(n)(Ims -∆Ψ(n)) -1 .
The eigenvalues of the second term are in the range [-∆ ψ, 0], and the eigenvalues of

Γ(n) are in the range of [1 -∆ ψ, 1].
Corollary 1: If f φ = 0 and K f ξ < 0, then there exists a ∆ > 0 for which P R is exact.

Proof: (Sketch) Observe that

k l=n Γ(l)
is equal to Ims plus terms that are norm-bounded by positive powers of ∆ ψ. Similarly,

τ k=n k l=n Γ(l)
is equal to (τ -n + 1)Ims plus terms that are norm-bounded by positive powers of ∆ ψ. We can make these terms arbitrarily small by choosing ∆ small. We therefore write

τ k=n k l=n Γ(l) ≈ (τ -n + 1)Ims.
Because f φ = 0, we have that

Φ(n) ≈ -(τ -n + 1)∆K f ξ .
This is strictly positive for all n ∈ N if K f ξ < 0.

Given ψ, which is easy to obtain from the system parameters, we can estimate the maximum value of ∆ that guarantees exactness. Given that ∆ is small, we need only ensure that Γ(n)K f ξ > 0 and not worry about terms with higher order products of Γ(n). This will be the case, for example, if we choose ∆ small enough that the smallest element of the vector K f ξ is greater than the largest magnitude element of ∆ ψK f ξ .

Corollary 2: If f ξ = 0 and f φ > 0, then there exists a ∆ > 0 for which P R is exact.

Proof: (Sketch) Following the same logic as Corollary 1, we can choose ∆ > 0 such that

Φ(n) ≈ Ims -∆K τ k=n J (ξ(k)) f φ . ( 11
)
The second term in the parentheses can be made arbitrarily small by choosing ∆ to be small. In this case

Φ(n) ≈ f φ , which is positive if f φ > 0.
A shortcoming of Corollaries 1 and 2 is that they can be limited to short time intervals. This is because if an interval is to remain constant, the number of time periods, τ , must increase as the step, ∆, decreases, and the approximations in the proofs of the corollaries do not hold for increasing τ . On the other hand, these are conservative sufficient conditions. For example, the second term in the parentheses of [START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF] will often be positive semidefinite or nearly so. This is because it is typical for K to be lower triangular with a negative diagonal and for J (ξ(n)) to be nonnegative and nearly diagonal. We therefore expect that exactness will sometimes hold for larger ∆ over longer time intervals.

We view these theoretical results not as a complete characterization of when exactness is guaranteed, but rather as evidence that P R is exact for a meaningful set of problems, and as guidance as to how to identify them. While there are certainly problems of interest for which their conditions do not hold, P R may nonetheless provide a useful and sometimes perfect approximation.

A. Steady state

It may be of interest to optimize (1) in steady state, e.g., when the solution does not change significantly on the timescale of interest, to find the best operating point, or to reduce the number of variables. We obtain a steady state optimization by dropping the time index and replacing the finite difference in (4c) with zero. The resulting (relaxed) optimization is:

P RS min ξ,ξ in ,T F(ξ, T ) (12a) such that T ≤ φ(ξ) (12b) 0 = V KT + N ξ + Ĉξ in (12c) ξ, ξ in , T ∈ Ω. (12d)
We remark that, in general, a solution to P RS is not guaranteed to be an equilibrium of (1). One special case in which guarantees do exist is the gradostat, which we discuss in Example V-A. We refer the reader to [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF] for a brief summary.

Theorem 3: P RS is exact if the network is outflow connected and at the optimal solution,

0 < Ims + K V N -1 J (ξ) -1 × K V N -1 ∇F ξ (ξ) -∇F φ (T ) . Proof: The Lagrangian of P RS is L = F(ξ, T ) + ρ (T -φ(ξ)) + λ V KT + N ξ + Ĉξ in .
Differentiating the Lagrangian and setting it to zero gives

0 = ∇F φ (T ) + ρ + K V λ J (ξ) ρ = ∇F ξ (ξ) + N λ.
We now solve for ρ. N is invertible due to outflow-connectedness. Then

λ = N -1 J (ξ) ρ -∇F ξ (ξ) , and
0 = ∇F φ (T ) + ρ + K V N -1 J (ξ) ρ -∇F ξ (ξ) .
Solving, we have

ρ = Ims + K V N -1 J (ξ) -1 × K V N -1 ∇F ξ (ξ) -∇F φ (T ) .
By complementary slackness, P RS is exact when ρ > 0.

As in the latter part of the previous section, we now assume that F ξ and F φ are linear with gradients f ξ ∈ R ms and f φ ∈ R rs .

Corollary 3: Suppose that for all ξ ≥ 0,

Ims + K V N -1 J (ξ) -1 ≥ 0, and
K V N -1 f ξ -f φ ≥ 0.
If at least one of the inequalities is strict, then P RS is exact. Like Theorem 3, Corollary 3 depends on the optimal solution, but to a lesser extent. Whereas Theorem 3 depends on ξ and T through the objective and J (ξ), Corollary 3 only depends on ξ through J (ξ). This is more manageable because J (ξ) is often positive and nearly diagonal. For example, if we assume assume that biomass is constant and all growth rates are Monod, as in Section III-B, then J (ξ) is positive on the diagonal and zero elsewhere.

V. EXAMPLES

A. The gradostat

The gradostat is a special case of (1) where in each tank i ∈ S, a single substrate, S i , is converted to a single type of biomass, X i . The conversion occurs at the rate φ i (S i , X i )/y, where y is the yield. Then ξ i = [S i , X i ] and κ i = [-1/y, 1] . In [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], several examples are given in which the gradostat is optimized over time and in steady state.

Here we first apply Corollary 1 to the gradostat. Suppose

F ξ (ξ(n)) = f S S(n) + f X X(n).
In this case, Corollary 1 holds if -f S /y + f X < 0. When dealing with the decontamination of undesirable substrates, we do not seek to minimize biomass, in which case f X = 0 and the condition is satisfied if f S > 0.

We now apply Theorem 3 to the gradostat in steady state. Let J S (ξ) ∈ R r×s be the Jacobian matrix of φ(ξ) with respect to S, and let J X (ξ) ∈ R r×s be the Jacobian matrix of φ(ξ) with respect to X. P RS is exact if

0 < Is + V N -1 -J S (ξ) /y + J X (ξ) -1 × K V N -1 ∇F ξ (ξ) -∇F φ (T ) .
It is straightforward to verify that when F ξ (ξ) = 0, this condition directly implies Theorem 1 in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF]. Similarly, we obtain Corollary 1 in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF] if we specialize Corollary 3 to the gradostat.

B. Wastewater treatment

In this example, we optimize an idealized wastewater treatment system consisting of three wastewater treatment plants of the city of Paris and its suburbs. The model is based on that in [START_REF] Robles-Rodriguez | A simple model of wastewater treatment plants for managing the quality of the Seine River[END_REF], [START_REF] Robles-Rodriguez | Management of an integrated network of wastewater treatment plants for improving water quality in a river basin[END_REF], which is based on the Activated Sludge Model no. 1 [START_REF] Henze | Activated sludge models ASM1[END_REF], and the influent data from the Inf rain 2006 dataset of [START_REF] Alex | Benchmark simulation model no. 1 (BSM1)[END_REF]. We implemented the model using the parser CVX [START_REF] Grant | CVX: Matlab software for disciplined convex programming, version 2.1[END_REF] and the solver Gurobi [START_REF] Optimization | Gurobi optimizer reference manual[END_REF] on a personal computer from 2014 with a 1.4 GHz dual-core processor.

Sewer networks typically have several valves and pumps for controlling flow rates. Modeling these as actuators would unfortunately introduce bilinearities, which, while amenable to techniques like McCormick relaxations [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I-Convex underestimating problems[END_REF], are out of our current scope. Here we instead assume that the flow rates are constant and fixed amounts of BOD and NH + 4 are to be allocated over the treatment plants in each time period.

The flow rates are given by Q in 1 = 0.1 m 3 /s, Q in 2 = 0.4 m 3 /s, and Q in 3 = 0.2 m 3 /s. All tanks have volume 1000 m 3 . In each plant i ∈ S, the state vector

ξ i = ξ BOD i , ξ NH + 4 i , ξ NO - 2 i , ξ NO - 3 i ∈ R 4
consists of biochemical oxygen demand, ammonia nitrogen, nitrite, and nitrate. The biomass concentrations in each plant and time period are assumed to be exogenous parameters, and therefore not components of the process state. This is an admissible assumption in the sense that the dynamics of the biomass concentration is typically much slower than that of the other process components, and has a larger amplitude. As a result of this assumption, it is not necessary to specify biomass influent concentrations.

In each tank i ∈ S, the elements of the process state have

kinetics φ BOD i , φ NH + 4 i , φ NO - 2 i
, and φ

NO - 3 i
. All assume Monod growth rates with parameters given in Table I, which comes from Table 1 in [START_REF] Ayed | Real-time optimisation of a wastewater treatment network[END_REF]. Because the biomass is constant, the growth kinetics can be represented as SOC constraints in P R . The stochiometric matrix for

κ i =       -1 0 0 0 0 -1 0 0 0 1/y NH + 4 ,NO - 2 i -1 0 0 0 1/y NO - 2 ,NO - 3 i -1       .
We used the implicit Euler method, Dn[ξ(•)] = (ξ(n)-ξ(n-1))/∆, with the stepsize ∆ = 1, which corresponds to 15 minutes. There are τ = 1345 time periods, so that the total time is two weeks. The boundary condition is ξ(0) = ξ(τ ). This could represent periodic operation, or exogenous conditions that are similar from week to week.

The process state must satisfy ξ BOD (n) ≤ 60 mg/L for each i ∈ S and n ∈ N . Without these constraints, most of the substrate would be directed to Plant 1, which is more efficient than the others. This constraint could represent, for example, regulatory limits on the pollution released by the plants.

In each time period, fixed quantities of BOD and NH + 4 , Ξ BOD (n) and Ξ NH + 4 (n), must be allocated over the treatment plants. These quantities are based on the Inf rain 2006 dataset of [START_REF] Alex | Benchmark simulation model no. 1 (BSM1)[END_REF], which covers 1345 15-minute intervals. Ξ BOD (•) is set to S S (readily biodegradable substrate), and Ξ NH + 4 (•) to S NH (NH + 4 and NH 3 nitrogen) in [START_REF] Alex | Benchmark simulation model no. 1 (BSM1)[END_REF]. The allocation is represented by the linear constraints

Ξ BOD (n) = 3 i=1 Q in i (n)ξ BOD,in i (n) Ξ NH + 4 (n) = 3 i=1 Q in i (n)ξ NH + 4 ,in i (n), for each n ∈ N . The other influent concentrations are ξ NO - 2 ,in i (n) = 3 and ξ NO - 3 ,

in i

(n) = 10 for i ∈ S and n ∈ N . The plant biomass concentrations are set to Xi (n) = 100 1 + (-1) i sin(10πn/τ ) for i ∈ S and n ∈ N . Observe that this has larger magnitude and varies slower than the other influents.

For each i ∈ S and n ∈ N , the optimization variables are ξ i (n), ξ BOD,in i (n), and ξ NH4,in i (n). The objective is to minimize the untreated wastewater released by the plants, given by n∈N i∈S

Q out i η i ξ i (n),
where we note that Q out i = Q in i , and η i = [2, 2, 0.3, 0.1, 0] . This objective was chosen to satisfy Corollary 1. Note, however, that the result does not fully apply due to the boundary condition, ξ(0) = ξ(τ ), and the upper limit on the process state.

The convex relaxation P R contains 145272 variables (in standard form) and took 17 minutes to solve. The 18 solver iterations accounted for only four seconds, and the rest of the time was used for preprocessing. Despite not fully satisfying Corollary 1, the solution was exact in all time periods. (n) to hit their concentration limits in Plants 1 and 3. When this happens, the remainder is allocated to Plant 2, which has little biomass at the time, and hence cannot transform the substrates as efficiently.

VI. CONCLUSION

We have formulated a convex relaxation for optimizing a broad class of bioprocesses. We proved that the relaxation is exact under simple conditions, and implemented it on a wastewater treatment example with over one hundred thousand variables. We believe that a wide range of problems can be approached in this manner due to the generality of the model and the tractability of SOCP.

One direction of future work is dealing with inexactness. Two options are deriving general convex underestimators to limit the relaxation gap, as in [START_REF] Taylor | Second-order cone optimization of the gradostat[END_REF], and finding local minima of the non-relaxed problem. In particular, the concave-convex procedure [START_REF] Yuille | The concave-convex procedure[END_REF] is wellsuited to the nonconvexity encountered here and would entail solving a sequence of SOCPs. We also intend to incorporate new elements into the model such as biomass death and recirculation, and to apply the relaxation in other contexts such as enzymes, where Michaelis-Menten kinetics [START_REF] Michaelis | Die kinetik der invertinwirkung[END_REF] have the same form as Monod kinetics.

Figures 1 and 2 Fig. 1 .

 21 Fig. 1. ξ in (•) from hour 175 to 275. The units are mg/L.

  i

Fig. 2 .

 2 Fig. 2. ξ(•) from hour 175 to 275. The units are mg/L.

Funding is acknowledged from the Natural Sciences and Engineering Research Council of Canada and the French LabEx NUMEV (Project ANR-10 LABX-20), incorporated into the I-Site MUSE, which partially funded the sabbatical of J.A. Taylor at MISTEA lab.