Assessing the Versatility and Robustness of Pore Network Modeling to Simulate Redox Flow Battery Electrode Performance

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.

Large-scale, stationary energy storage technologies, such as redox flow batteries (RFBs), are attracting widespread interest due to their ability to compensate for the intermittent nature of renewable energy sources (i.e., solar photovoltaics and wind power). [1][2][3] RFBs are rechargeable electrochemical systems that interconvert chemical and electrical energy by leveraging redox couples, dissolved in liquid electrolytes, that are pumped through an electrochemical stack to charge and discharge the battery. 4 Their most appealing features are the ability to decouple power and energy capacity which is not possible with sealed batteries (e.g., lithium-ion), high round-trip efficiency, extended durability, fast response times, low environmental impact, and geographic independence. 5,6 Despite their promise, RFBs have so far seen limited market penetration due to their current elevated costs and suboptimal performance, motivating substantial research efforts on alternative redox chemistries and novel reactor architectures. [7][8][9][10][11][12] As a core component of the flow battery stack, the porous electrode largely determines the battery performance by providing both the active surface area for electrochemical reactions, the solid architecture for electron conduction, and the open porosity for electrolyte flow and mass transport. Hence, optimizing the electrode microstructure and surface chemistry is an effective approach to decrease reactor and system costs.
Electrode engineering has been driven primarily by empirical approaches, targeting increases in available surface area, kinetic activity, and wettability for the most studied systems (e.g., all-vanadium), 7,13,14 but this approach is time-and resource-consuming. Furthermore, each redox chemistry and reactor concept necessitates a unique set of electrode properties, which results in a very broad design space. Recent advances in computational sciences enable the simulation of complex electrochemical flows within porous electrodes, which can be leveraged to understand the role of the porous electrodes and ultimately to accelerate their design pipeline. In this context, simulations have been deployed for RFBs at multiple length scales. 15,16 For example, cell-level macroscopic continuum models have been employed, [17][18][19][20][21][22][23][24] aided in investigations of cell-averaged properties, 18 proved to be effective for predicting the cell performance for multiple charge-discharge cycles for all cell designs. 21 However, these models build on volume-averaged grids and fail to capture the complex effects of the electrode threedimensional structure (e.g., local variations in porosity or surface area, inhomogeneous pore size distributions) which results in inaccurate predictions of the reactor performance. Therefore, to accurately describe the three-dimensional structure of the material, microstructure-informed simulations have been deployed.
Highly detailed models, often referred to as direct numerical simulations, using finite element and lattice Boltzmann methods have been employed to gain insight into the transport within porous materials. [25][26][27] These models solve the velocity field of the electrolyte using the voxels of a tomogram as the mesh, and can therefore deal with highly complex geometries, such as carbon fiber-based electrodes. Qiu et al. employed a multiphysics lattice Boltzmann method to obtain the local velocity field at the pore-scale 26 and investigated the effect of the electrode porosity and electrolyte flow rate on overall RFB performance, quantifying the effect of fuel starvation within vanadium systems. 27 The authors found that low porosity electrode structures result in a more uniform current density and lower absolute overpotential distributions at the expense of increased pumping requirements. Recently, Bao et al. coupled a multi-scale model with a deep neural network to optimize a timevarying electrolyte inlet velocity for an all-vanadium RFB leading to a 74% reduction in the pumping power consumption. 28 Additionally, Yoon et al. investigated the effect of the electrode local porosity on battery efficiency, showing that lower porosity at the electrolyte inlet regions resulted in superior energy efficiency as it provides a larger active area where the reactant concentration is highest. 29 These studies demonstrate the potential of three-dimensional modeling to inform electrode and reactor design; however, they require large computational power to solve the defined set of equations. 30 The z E-mail: a.forner.cuenca@tue.nl = Co-first authors. *Electrochemical Society Member. large number of nodes required to model even a few pores and the non-linearity of the physics limits the application of these models to a small subdomain of the electrode. Thus, to accurately describe reactor performance at larger sizes, modeling approaches that can simultaneously capture electrode microstructures and simulate varying phenomena at the reactor level need to be developed.
In recent years, pore network modeling (PNM) has been employed as an alternative modeling technique that can capture microstructural effects at the mesoscale at a reasonable computational cost by invoking geometrical simplifications on the porous network. 9,10,31 In PNM, the complex porous space is approximated by a network of spheres (pores) connected by cylinders (throats) and one-dimensional analytical solutions are leveraged to compute the transport of mass and charge within the network. Agnaou et al. demonstrated that this simplification of the model equations results in a computational time reduction up to 10 4 times compared to finite element simulations within the same geometry, with a limited loss in computational accuracy. 30 PNM has been employed in the simulation of multiphase flows of fuel cell gas diffusion layers [32][33][34] and has been recently applied to porous electrodes in RFBs. Sadeghi et al. built a PNM framework to study the liquid side of the hydrogenbromine flow battery. 10 The authors investigated the effect of the electrode porosity on the performance of a simulated SGL 25AA electrode and computed a relationship between electrode permeability and voltage efficiency. 10 Lombardo et al. presented a framework for transient PNM to simulate the coupled flow, species, and charge transport within a vanadium flow battery using an extracted Toray 090 electrode. 9 They displayed the existence of masstransport limited areas that significantly contribute to performance losses. Finally, other groups have employed PNM as a tool to compute the permeabilities and parasitic losses in carbon papers and cloths from tomographic reconstructions of dry electrodes. 35,36 These previous studies offer an optimistic prospect for the use of PNM as a platform for electrochemical simulations of porous electrodes in RFBs. However, there are notable knowledge gaps that must be addressed to deploy PNM broadly. First, methods that can reliably represent complex electrode geometries while affording larger sample sizes (i.e., cell level) are lacking as most research has focused on a small domain of the electrode. 9,10,35,36 Second, the investigation of geometrical algorithms that afford large versatility and robustness across microstructures and redox chemistries has not been performed. Inspired by these knowledge gaps and previous work in the field, we seek to answer the following questions: "can we build a framework that is robust and versatile to accurately describe different electrode structures and electrolytes?" and "what is the minimum number of fitting parameters required to accurately describe the cell performance with PNM?" Furthermore, we aim to leverage this framework to theoretically investigate and understand the performance differences between cloths and papers under different electrolyte systems, which has recently triggered significant attention of experimentalists. 7,[37][38][39] Here, we introduce a computationally inexpensive, chemistryagnostic, and microstructure-informed PNM framework using the open-source software OpenPNM. 31 An iterative algorithm enables the coupling of both half-cells and utilizes a network-in-series approach to account for species depletion over the length of the electrode. We study two fibrous electrodes -a Freudenberg paper and an ELAT Cloth -and two electrolytes -an Fe 2+ /Fe 3+ couple in aqueous media and a TEMPO/TEMPO + couple in a non-aqueous media -to screen a broad property space and, in turn, assess the versatility of the modeling framework. First, X-ray tomographic microscopy is performed to obtain dry microstructures. The network extraction is evaluated by comparing the pore size distributions of the extracted pore networks with mercury intrusion porosimetry. Second, the model is validated against experiments, using single-electrolyte flow cells to enable precise quantification of the cell performance for both electrodes with the non-aqueous electrolyte. Third, we compare simulations to flow cell experiments using an aqueous electrolyte and assess the versatility and robustness of the model revealing that incomplete wetting of the porous structure impacts physical descriptors used for single-phase flows. Fourth, a systematic analysis of the electrochemical performance in tandem with local properties (i.e., pressure drop, current density, reactant depletion, and overpotentials) is performed to elucidate microstructure-performance relationships. Finally, a parametric study was performed to identify efficient operation envelopes.

Model Development
The present study offers a numerical platform that can be leveraged to simulate the local transport within porous electrodes for different RFB chemistries at a low computational cost. The modeled domain (Fig. 1a) consists of two symmetric, mirrored 1 mm long by 1 mm wide porous electrodes, extracted using X-ray tomographic microscopy (XTM). Symmetric flow cell simulations were conducted with co-current operation of the anodic and cathodic half-cells. In this setup, the oxidation reaction of the redox couple occurs in the anodic half-cell and the reduction reaction of the same redox couple takes place in the cathodic half-cell. Local transport equations are solved within the porous electrodes, while the location of the flow channels, current collector, and membrane are dictated by the boundary conditions. Coupling of the two half-cells is achieved using the boundary conditions of the membrane, which is treated as a macro continuum entity and thus only the overall macroscopic ionic resistance of the membrane is considered.
Pore networks (PNs) are an idealization of real tomographic extracted porous structures. In PNMs, it is assumed that the void space within the porous structure can be approximated by spherical pores and cylindrical throats. The bulk solution in each pore is assumed to be well-mixed within the pore body and transport in the PN occurs via the throats. This idealization of the void space allows for the reduction of complexity of the considered model equations, while still retaining information about the microporous structure and topology of the porous electrode. 31 Consequently, the computational time to run one network with the network-in-series approach, for applied potentials of −1 to 0 V (to simulate the consumption of reactant species, we run the model in discharge mode) using −0.05 V step intervals, was only ∼100 min for a 1.7 cm long electrode using a single Intel® Core TM i7-8750H CPU, depending on the electrolyte chemistry and velocity.
Model equations.-The considered physics within the threedimensional porous electrode and electrolyte account for fluid transport, species transport, reaction kinetics, and conservation of charge. The flow fields, electrodes, and membrane are coupled by the applied boundary conditions. The governing equations and the formulated boundary conditions within the electrochemical cell are described for the case of cell discharge and are based on the convention of describing the flux of positive charges. To simplify the complexity of the model, the following assumptions are made: 1. The electrochemical cell operates at steady-state, iso-thermal conditions. 2. The dilute-solution theory is applied and only the active species are included. 3. The crossover of active species through the membrane is neglected. 4. The electric potential of the solid electrode phase is assumed to be constant due to its high electronic conductivity. 32 5. The flow within the porous structure is classified as a creeping flow, due to the small flow dimensions and low Reynolds numbers, thus inertial effects were neglected when computing the pressure drop. 6. The charge flux within the liquid phase is approximated using the bulk conductivity of the electrolyte, neglecting migration effects. 7. The half-cell reaction is described as a single-step reaction mechanism and neglects the formation of intermediates that are adsorbed on the electrode surface and side reactions. 8. The electrochemically active surface area of a pore is equal to the extracted area obtained with the SNOW algorithm, which utilizes the marching cubes algorithm to obtain areas of smoothed surfaces from the voxelated image. 9. The electrochemical reaction is only occurring in the pores and not throughout the connecting throats, though the area of the throats was assigned to the pores during the network extraction.
The fluid transport within the porous electrode is described by the steady-state Navier-Stokes equations. The total mass balance for the electrolyte phase around pore i in the network is defined as where ρ is the electrolyte density, u ij the fluid velocity from pore i to pore j, S ij the cross-sectional area of the throat connecting pore i and j, and N T the number of throats connected to pore i. By treating the throats as cylindrical tubes, the Hagen-Poiseuille equation can be used to express the velocity field based on the pore pressures as described by where p is the pressure in pore i or j, and α ij is the hydraulic conductance of the throat connecting pore i and j where μ is the dynamic viscosity of the electrolyte, and L ij the conduit length of the connecting throat.
The mass transport of chemical species in the porous electrode is described by the advection-diffusion-reaction equation. The conservation equation around each pore for the active species that is consumed is described by where m ij is the mole flux of the active species between pore i and j which contains both advective and diffusive terms, R i the net reaction rate of the active species, λ i the dimensionless stoichiometric coefficient in the corresponding half-reactions, I i the applied current in pore i, n the number of electrons participating in the halfreaction, and F the Faradaic constant. The mole flux is derived from the exact solution for the one-dimensional advection-diffusion equation, based on the exponential approach described by Sadeghi et al., 40 defined as where c is the concentration of the reacting species in pore i or j, and = Pe u L D ij ij the local Peclet number at the considered throat, with D the diffusion coefficient of the active species in the electrolyte.
The applied current in both the anodic and cathodic compartments is described by the Butler-Volmer equation defined as where I i is the applied current in the anodic or cathodic compartment, j 0 the exchange current density in the anodic or cathodic compartment, A i the electrochemically active internal surface area of pore i, c i s , the concentration of the reduced or oxidized form at the electrode surface, c ref the reference concentration of the reduced or oxidized form at which the exchange current density was measured, η the overpotential in the cathodic or anodic compartment, R the universal gas constant, T the operating temperature, and α the anodic or cathodic reaction transfer coefficient, where α α + = n a c with n the number of electrons transferred in the electrochemical reaction (for all tested chemistries in this study = n 1). The overpotentials in both half cells are given by a s e where φ s is the solid phase potential, φ e the liquid phase potential, and E oc the equilibrium potential defined against the same reference electrode. The concentration of the active species at the electrode surface is estimated from the corresponding bulk concentrations through the mass transfer of the active species towards the electrode surface where k m is the mass transfer coefficient of the considered species, and c i b , the bulk concentration of the active species. Combining Eqs. 6 and 10, and 7 and 11 by substituting for the surface concentrations, the Butler-Volmer equation for a one-electron electrochemical reaction including mass transfer limitations is obtained. It may be verified that the resulting equations reduce to the case without mass transfer limitations (Eqs. 10 and 11) when the mass transfer terms approach infinity ( → ∞ k m ). An estimation of the mass transfer coefficient towards the porous electrode was obtained by neglecting inertia effects and applying the film theory, assuming that the film layer is equal to the pore radius of every pore where d p is the pore diameter. The conservation of charge around each pore is coupled to the species transport in the electrochemical cell by the current and the potential field described by where I i is the applied current defined by the Butler-Volmer equation, and I ij the charge flux from pore i to j, which is proportional to the potential field in the liquid phase ij is the electrical conductance of the connecting throat, σ l the bulk electrolyte conductivity, and φ the liquid potential in pore i or j.
Boundary conditions.-For the fluid transport boundary conditions (Fig. 1b), the inlet pressure boundary condition is determined by setting a target inlet velocity. From the target inlet velocity, the target flow rate is calculated by where Q in is the flow rate of the electrolyte, u in the inlet velocity of the electrolyte, and A in the geometrical inlet area of the electrode. Subsequently, the inlet pressure at the bottom boundary pores is computed such that the total flow rate entering the network matches the target flow rate. The discharge pressure was set to zero. Furthermore, no-flux boundary conditions were imposed at the width and thickness of the electrode x x x x x and z 0 , , , , At the current collector-electrode interface in the cathodic compartment, the boundary condition for the solid potential is equal to the given cell voltage, while in the anodic compartment the solid potential is equal to zero as a result of the symmetry of the modeled domain where V cell is the given cell potential. The boundary condition for the electrolyte potential at the membrane in a specific half-cell is iteratively calculated from the electrolyte potential at the membrane in the other half-cell, using Ohm's law to include the average voltage loss across the membrane interface where φ Δ m is the average voltage loss across the membrane interface, R m the resistance of the membrane, I m the current passing the membrane, and φ m the electrolyte potential at the membrane in the anodic or cathodic compartment. At the other boundaries, the boundary condition of the voltage was defined as Iterative algorithm.-The utilized numerical algorithm consists of the coupled mass and charge transport within the porous electrode. In the developed algorithm, the individual transport equations within the anodic and cathodic compartments are updated sequentially within an iterative algorithm. The detailed flowchart of this algorithm is illustrated in Fig. 2.
Due to the assumption of a dilute electrolyte, the viscosity and density of the fluid are assumed constant so fluid transport within the porous electrode is considered to be independent of the other transport processes. Therefore, it can be directly solved to obtain the velocity and pressure fields with Eqs. 1 and 2. The velocity field of the liquid electrolyte is provided to the iterative algorithm to solve the species and charge transport equations, which consist of two nonlinear systems of equations (Eqs. 4 and 13) that are coupled by the Butler-Volmer equation (Eqs. 6 and 7). The iterative scheme consecutively solves the advection-diffusion-reaction equation (Eq. 4) and the potential field (Eq. 14). Both half cells are coupled, using the found solution for the electrolyte potential at the membrane in one-half cell to calculate the electrolyte potential at the other side of the membrane in the other half cell. In this step, Ohm's law was used to account for the membrane resistance (Eq. 21). The calculated electrolyte potential is subsequently used as a boundary condition for the charge transport equation in the second half cell.
The numerical model solves the coupled transport equations over a given voltage range starting from the open-circuit voltage (0 V for the performed simulations), with an initial guess for the overpotential in the first iteration of 0 V. The initial concentration within all pores was set to the inlet concentration of the active species. While progressing to higher absolute cell potentials, the initial guesses for the concentration, potential field, and overpotential are based on the solution of the previous cell potential.
Numerical convergence is achieved when the specified relative and absolute tolerances are met. These tolerances are based on the total current generated or consumed in a half-cell defined as where I tot is the total current generated or consumed in a half-cell, I loc the (local) current generated in each pore, and N p the number of pores. The relative error can be calculated for both compartments by where ϵ rel is the relative error, and n the iteration number, with a maximum of 5 × 10 4 iterations. The absolute error, ϵ , abs consists of the difference in total current between both half-cells and is defined as The tolerances for the relative and absolute errors were set to 5 × 10 −5 and 6 × 10 −4 A cm −2 , respectively. To counteract the divergence of the solution due to the highly nonlinear nature of the system, two numerical strategies were employed. The first strategy consists of under-relaxation during the updating of the concentration and potential fields. 10 In under-relaxation, a relaxation factor is used to dampen the obtained solution where X is the concentration or the potential of species i, and ω is the relaxation factor, set to 0.1. The second strategy consists of linearization of the charge transport source term. When the considered active species has facile kinetics, the mass and charge transport equations become strongly coupled. In this case, linearization of the current generation term can mitigate divergence. The derivation of the source term is further elaborated on in the Supporting Information (Section S1 (available online at stacks.iop. org/JES/169/040505/mmedia)).
Network-in-series model.-To account for species depletion, the numerical framework was modified to match the entire length of the laboratory electrode (1.7 cm) using a network-in-series approach. In this approach, multiple PNs are placed behind one another, where the concentration at the end of the previous network is considered as the inlet concentration of the next PN. To map the concentration between the boundary pores of the two networks, each network was manipulated to become a mirrored copy (mirrored in the flow dimension) of the previous network (see Fig. S1). The mass and charge transport algorithms are solved independently for each anode-cathode couple in the series.

Experimental
Electrode materials.-Two different commercial carbon fiberbased electrodes were investigated in this work. A carbon paper, Freudenberg H23 (Fuel Cell Store), and a woven carbon cloth, ELAT Cloth electrode (Fuel Cell Store). Freudenberg H23 has an uncompressed thickness of 210 μm and a porosity of 80%, and the ELAT Cloth has an uncompressed thickness of 380 μm and a porosity of 82%. 7 Electrode density measurements were performed for six different electrode pieces for both electrode types (1.7 cm × 1.5 cm × 210 μm or 375 μm) using a Sartorius ED224S analytical balance (± 0.1 mg uncertainty). The electrodes were used as received. In addition, three ELAT Cloth electrodes were functionalized with electrografted taurine (2-aminoethane-1-sulfonic acid) to investigate the impact of electrode hydrophilization on the performance and model validation.
Flow cell experiments.-Symmetric flow cell experiments were conducted in a laboratory-scale flow cell platform for both the aqueous and non-aqueous electrolytes (see Fig. 1a for a schematic representation of the cell). The corresponding redox reactions are 5 where E 0 is the standard reduction potential. Separate tubing was used to connect the single solution reservoir to the cell in-and outlets to ensure a constant state of charge at the cell inlets. Graphite current collectors milled with a flow-through flow field were employed to force convective electrolyte flow through the porous electrode, resulting in near unidirectional flow parallel to the membrane. This simplified design is well suited to perform fundamental investigations on the role of the electrode microstructure and electrolyte properties without the complexities of spatially varying flow fields (e.g., interdigitated, serpentine). One electrode with an external area of 2.55 cm 2 (1.7 cm × 1.5 cm) was placed on top of each current collector using incompressible polytetrafluoroethylene gaskets with a total thickness of 210 μm or 375 μm, respectively for the Freudenberg H23 and the ELAT Cloth, ensuring minimal compression of only a few microns. Here, unless otherwise stated, we elected to study flow cell performance relationships without compressive effects. Thus, both the tomograms and the electrochemical tests were performed using the same, uncompressed, conditions to enable representative comparisons. We acknowledge that practical systems employ compressive forces (20%-50%) to minimize contact resistances, which are accounted and corrected for in our study, and that the electrode microstructure changes as a result of compression. 42 However, dedicated investigations on compressive strengths are beyond the scope of this study and have been recently treated elsewhere. 43,44 Given the small area of the cell and relative thickness of the endplates, this minimal compression was sufficient to maintain reliable operation. A cation-exchange Nafion 211 membrane (Fuel Cell Store, 25.4 μm), was used for the Fe 2+ /Fe 3+ system, and a porous separator, Daramic 175 (SLI Flatsheet Membrane, 175 μm), was used for the TEMPO·/TEMPO + system. Subsequently, the cell was tightened with a torque-controlled screwdriver to 2.2 N m. To account for differences in electrode thickness, the performance was compared at the same inlet velocities of 1.5, 5.0, and 20 cm s −1 , calculated by dividing the volumetric flow rate set by the pump by the geometrical inlet area of the flow-through flow field (thickness multiplied by the width of the electrode). The electrolyte was pumped through the cell using a Masterflex L/S® Easy-Load® II pump and LS-14 tubing in co-current operation of the half cells to promote the wetting of the porous carbon electrodes. Electrochemical experiments were conducted in charge mode, performed inside a nitrogen-filled glove box for the non-aqueous electrolyte, and outside the glove box for the aqueous electrolyte. Polarization curves were obtained by employing constant voltage steps of 0.05 V for 1 min and measuring the steady-state current for a lead-corrected voltage range of 0 to 1 V. For the Fe 2+ /Fe 3+ system an Ivium (Ivium-n-Stat) potentiostat was used and for the TEMPO/TEMPO + system a Biologic VMP-300 potentiostat. The polarization curves were repeated three times for new assemblies and new electrolyte solutions for each flow rate in descending flow rate order.
Electrolyte density, viscosity, and conductivity measurements.-Electrolyte density measurements were performed in triplo for the two electrolytes using a Sartorius ED224S analytical balance and an electrolyte volume of 10 mL. The electrolyte viscosity was assumed equal to the viscosity of the solvent. Conductivity measurements were performed using a two-electrode custom conductivity cell, similar to the setup used in Milshtein et al. 45 The compartment of the conductivity cell was flooded with electrolyte and sealed shut. Electrochemical impedance spectroscopy (EIS) was conducted using the Biologic VMP-300 potentiostat at open-circuit voltage and room temperature (20°C) with an amplitude of 10 mV and a frequency range of 1 kHz−100 kHz, 8 points per decade, 6 measurements per frequency, and a waiting time of 0.10 period before each frequency. The high-frequency intercept was identified as the value of the total resistance (cell, lead, and electrolyte resistance). EIS measurements were performed in triplo and were repeated five times for new assemblies and new electrolytes. A calibration curve was obtained using aqueous conductivity standards (0.01 M, 0.1 M, and 1.0 M aqueous potassium chloride) and an empty cell was measured to obtain the combined cell and lead-resistances.
Resistance measurements.-EIS measurements were performed at open-circuit voltage where the high-frequency intercept was identified as the value of the flow cell resistance capturing e.g., the membrane resistance and electronic contact resistances. The flow cell resistance (high-frequency intercept corrected for the contact and lead-resistance) for the Fe 2+ /Fe 3+ chemistry with the Nafion 211 membrane and NaCl supporting electrolyte was determined to be 0.39 ± 0.02 Ω, while it was 0.77 ± 0.01 Ω for the TEMPO . /TEMPO + chemistry with the Daramic 175 membrane and TBAPF 6 electrolyte. The obtained values were used in the model as membrane resistivity values, by multiplying the high-frequency intercept resistance by the external area of the electrode (2.55 cm 2 ).
Exchange current density.-The exchange current density ( j 0 ) was extracted from the experimental data by fitting the low current density region (< 100 mA cm −2 ) of the polarization curve to the Butler-Volmer Eqs. 6 or 7): where i is the measured current density, and η act the activation overpotential obtained after iR Ω -correction of the experimental data and assuming mass transfer overpotentials to be negligible, resulting in c i s , being equal to c . ref The exchange current densities were obtained for each electrolyte at all inlet velocities and for both electrodes.
X-ray tomography and image processing.-Electrode materials were scanned in an uncompressed state using a laboratory micro-CT (Scanco Medical μCT 100 cabinet microCT scanner, holder type U50822 with a diameter of 9 mm and a height of 78 mm) at an isotropic resolution of 3.3 μm/voxel. The scans were carried out using a peak potential of 55 kVp, a current of 72 μA, 4 W, and a 0.1 mm aluminum filter to acquire 312 projection images over 360 degrees.
Gray-scale images produced using XTM contain a certain degree of noise that has to be reduced. Care was taken to preserve phase edges while applying noise reduction filters to limit losses of features in the virtual model. In this work, a two-dimensional median filter with a radius of 2.0 pixels was used to reduce the noise in the grayscaled image. Subsequently, each voxel in the gray-scale image has been assigned to either the solid or void phase using a K-means cluster segmentation filter. 46 The reported median filter and K-means cluster segmentation filter were imported as plug-ins in the image processing software ImageJ. 47,48 To prevent boundary effects from sample cutting, a 1.0 mm high and wide sample selection was selected from the center of the processed image, where the thickness of the electrode was reduced by the removal of the first and last few image layers in this dimension to prevent artifacts during the pore network generation.
The electrode geometry used in the PNM was generated by the extraction of the physical electrode structure from XTM images, in combination with image processing steps consisting of smoothing, segmenting, and network extraction by the SNOW algorithm, using the equivalent diameter. 49 This open-source algorithm was developed for the extraction of porous networks from tomographic images and was validated for porous media with a wide range of porosities (20%-85%). 49 The network extraction operation was performed using a single Intel® Core™ i7-8750H CPU. For the extracted networks, a standard deviation value of the applied convolution mask of 0.4 and a value for the radius of the spherical structure element of 5 voxels were used. 49

Results and Discussion
Network extraction.-The reliability of the modeling framework largely depends on the ability of the simplified porous network geometries to represent the complex microstructure of fibrous electrodes. Thus, before validating the simulated electrochemical performance, we first assess the accuracy of the network extraction procedure by comparing relevant microstructural parameters to conventional methods and literature values.
The extracted networks are represented graphically in Figs. 3a -3b, including an overlay of the spherical pores over the original tomograms. As a result of the organized fiber bundles of the woven electrode, a bimodal pore size distribution (PSD) is observed consisting of small pores in between the fibers and larger pores at the intersection of the fiber bundles (Fig. 3b). On the contrary, the non-woven material contains a unimodal distribution. To quantify this, we compute the PSDs of the networks and compare them with mercury intrusion porosimetry (MIP) data (Figs. 3c-3d) obtained in a previous study. 7 The exact steps taken to evaluate the intrusion data can be read in the Supporting Information (Section S2) and it should be noted that MIP relies on intruding mercury onto the electrode pores and assuming a perfect cylindrical shape. Hence, the resulting comparisons are qualitative but give a good indication of the quality of the extraction procedure.
The pore size distributions represented in Figs. 3c-3d reveal a stark microstructural difference. A good agreement was observed for the PSD of the paper electrode obtained by both methods (Fig. 3c), which features an average pore diameter of ∼15 μm. However, pores larger than 60 μm are absent in the extracted network as a result of both the presence of artificial bigger pores between stacked electrodes in the MIP measurements, and the required image trimming steps taken in the extraction procedure (the surface layer of conventional fibrous electrodes is characterized by a higher porosity and a lower fiber density packing than the bulk electrode 49 ). The extracted PSD of the cloth electrode, on the other hand, differs from the MIP data (Fig. 3d). First of all, the distribution is shifted towards larger pore sizes, as the void space between the individual continuous fibers in the threads of the cloth cannot be accurately represented by pores and throats. Secondly, similar to the paper electrode, the SNOW extraction of the cloth results in the construction of a network with a maximum pore size of 260 μm for the larger pore segments, while the intrusion data displays bigger pores. To investigate whether inhomogeneities in the larger pore segment of the extracted network also play a role, a bigger network (170% of the original size) was extracted (Fig. S2). This extracted network matches the larger pores observed with MIP better. Nonetheless, an investigation into the performance of this bigger extracted network showed only a slightly better match with the experimental data, but similar curvature offsets were observed for the smaller extracted network (Fig. S3, Fig. 4b, and 5b).
The geometrical surface area, bulk electrode porosity, in-plane and through-plane permeability, anisotropy ratio, and thickness of the extracted PNs are compared to literature values in Table I. The in-plane permeability is defined as the average permeability in the direction parallel to the membrane (x and y-directions), while the through-plane permeability is defined as the permeability perpendicular to the membrane (z-direction). The anisotropy ratio is defined as the ratio between the in-plane and through-plane permeabilities and gives information about the pore alignment. A slightly reduced permeability, thickness, and porosity were observed for all the extracted networks compared to the experimental values. A larger deviation is observed between the extracted networks' internal surface area and the literature obtained value, with an increase in the surface area of 2.2 × for the paper network and of 9.5 × for the cloth network compared to the literature. The explanation of the higher internal surface area for the extracted networks is a result of assumptions taken in the experimental determination of the electrochemically active surface area (ECSA) through analysis of the electrochemical double-layer capacitance, such as the assumed specific conductance of the electrode, as they can lead to underestimations in the experimentally obtained ECSA. 8, 50 Tenny et al. investigated four volume-specific surface area techniques (MIP, Brunauer-Emmett-Teller, electrochemical double-layer capacitance, and morphological methods) for woven electrodes and found that the ECSA value estimated from the electrochemical double-layer capacitance (0.05-0.08 × 10 5 m 2 m −3 ) was several magnitudes lower than estimated by the other methods (1.1-2.6 × 10 5 m 2 m −3 ), which are in-line with our obtained pore network ECSA. 37 Therefore, the surface area was additionally obtained from the tomograms using the open-source image analysis package PoreSpy, which only relies on the geometrical surface area instead of the ECSA, and the surface roughness is only partially captured as a result of the XTM systems resolution. Nevertheless, the extracted networks differ by a factor of 1.5 with this estimation of the surface area as a result of the translation of the void space to pores and throats in the pore network. Despite these intrinsic challenges with surface area determinations, we remain confident about the extracted values since we determine the kinetics of the electrode-electrolyte couples from the experimental data in terms of the exchange current density, which corrects for any errors in the surface area.  Model validation.-We then assess the accuracy of the electrochemical model with experimental data for the non-aqueous electrolyte (the electrolyte properties are listed in Table II) and both electrode microstructures, by evaluating the single electrolyte cell polarization curves (Fig. 4). Importantly, the data presented in Fig. 4 only uses the obtained exchange current density extracted from the experimental measurements (Table SII) to obtain the electrolyte-electrode kinetics as an electrode-dependent parameter. Kinetic information needs to be provided to the PNM as the model is focused on the mesoscale, missing the capability to extract kinetic phenomena taking place at the molecular scale. Therefore, we elect to extract the kinetic information from experimental data, which is a common practice in electrochemical engineering modeling. 55,56 The exchange current densities obtained from the experimental data range from 1870-460 A m −2 for the non-aqueous chemistry and are in the same range as reported in the literature, using other electrode structures (375 A m −2 for the TEMPO · /TEMPO + chemistry 57 ). The extracted exchange current densities are depending on the electrolyte flow rate as within the timeframe of the experiments, it is likely that (fully wetted) steady-state conditions were not reached and that therefore the accessible surface area varies with different flow rates, resulting in velocity-dependent exchange current densities.
The PNM framework can predict the electrochemical response of a kinetically facile fully wetted system (non-aqueous) without the use of model-fitting parameters. This can be observed as the comparison between the model and experimental response for the non-aqueous electrolyte for both electrode structures, shown in Fig. 4, together with the root-mean-square error and the mean relative error (Table SIII), revealed a good match. The extracted exchange current densities provide a more accurate representation of the kinetics of the specific electrode-electrolyte systems; however, when running the simulations with the literature value (375 A m −2 ), the model was still able to capture the experimental response (Fig. S4). Thus, the proposed modeling framework was validated for the fully wetted non-aqueous electrolyte system with a symmetric cell design.
Incomplete wetting and parameter adaptation.-After the model validation with the non-aqueous electrolyte, the versatility and robustness is evaluated by comparing the simulated performance with experimental data for an aqueous electrolyte (the electrolyte properties are listed in Table II). The exchange current densities obtained from the experimental data range from 652-98.7 A m −2 for the aqueous chemistry and are again in the same range as reported in the literature (23.0 A m −2 for the Fe 2+ /Fe 3+ chemistry 58 ).
Whereas the behavior of the non-aqueous electrolyte is wellcaptured by the model, the simulated performance of the aqueous electrolyte shows significant deviations from the experimental data (Figs. 5a and 5b, Table SIII) which is hypothesized to be a result of incomplete wetting of the internal porous structure in the experimental measurements. The developed modeling framework assumes that the electrode is fully wetted, meaning that the entire electrode's geometric surface area is electrochemically active and that fluid transport in the open space of the porous structure can be characterized as single-phase flow. However, it has been shown that during the operation of RFBs, unfavorable interactions between aqueous electrolytes and carbon electrodes can result in gas hold-up within the porous structure. 8,[59][60][61][62] Especially at lower velocities this may lead to liquid channeling, which in turn leads to incomplete electrode utilization and hence decreased cell performance compared to the simulated performance assuming fully wetted electrodes. Due to their lower surface tension, organic electrolytes (i.e., acetonitrilebased electrolytes) have been demonstrated to fully wet carbon fibrous electrodes. 7 It is hypothesized that for this reason, the experimental performance of the TEMPO · /TEMPO + electrolyte is well-captured by the model, in contrast to the aqueous electrolyte.
The wetting behavior within RFB operation involves a complex interplay of concentration depletion, surface area utilization, and preferential pathway channeling. However, in this framework, the impact of these phenomena is not individually accounted for but encapsulated within the mass transfer coefficient, k . m The inability to capture the effect of these individual phenomena is an important limitation of this PNM compared to higher detailed computation models, such as lattice Boltzmann methods. 62,63 To account for the sub-optimal wetting in aqueous systems, we propose a fitting procedure for k m using SciPy's least squares method to minimize the error between the data points obtained from the experimental dataset and the model. 64 The electrochemical response of the corrected paper electrode system (Fig. 5c) yields a good agreement with the experimental data. It is observed that a greater k m correction factor was required (diverting from unity, factors fluctuate between 0.11-1, Table SIV) when operating with a lower hydraulic driving force (i.e., lower electrolyte velocity), as higher electrolyte velocities favor greater pressure drops, and thus a higher electrode saturation. This result supports the claim that the aqueous electrolyte, with a contact angle around 80°−90°with carbon fiber electrodes, 8 incompletely wets the carbon fibrous electrode in the experimental measurements leading to limited utilization of the available electrode geometric surface area. The corrected k m value for the cloth electrode however does not result in as good of a fit between the simulated and experimental data (Fig. S6), although better than the uncorrected data. The cloth electrode is more likely to suffer from contact resistances at minimal compressions compared to the flat paper electrode, as the inhomogeneous weave pattern of the cloth reduces the contact between the current collector and the electrode. At normal cell assembly conditions, the electrode is compressed to a certain degree, generally 20-30% 65 to minimize contact resistances, without compressing the electrode microstructure too much to prevent an increased pressure drop. 44 Therefore, impedance measurements were conducted to evaluate the medium frequency ranged arc, roughly translated to the charge transfer overpotential in the electrochemical cell. 18 The medium frequency arc of the uncompressed electrode displays a significant velocity dependency (Fig. S7a). This indicates that, in addition to small changes in activation overpotential with electrolyte velocity affecting the concentration of active species, the geometric surface area is utilized to different degrees. We hypothesize that, under minimal compression, the electrolyte could flow by the woven electrode and that this is affected by the electrolyte velocity, inducing different pressure drops and electrode utilization. To further investigate this hypothesis, we perform measurements with a 48% compressed woven electrode (Fig. S7b) and find that the medium frequency arc does not show a significant velocity dependency as the electrolyte is now more likely to be forced through the electrode. Nevertheless, both uncompressed (∼ 0.8-1.4 Ω for all tested velocities) and compressed (∼1.1-1.3 Ω for all tested velocities) cells feature a pronounced medium frequency arc in the EIS data as a result of the incomplete wetting causing partial access of the geometrical surface area. To further study the influence of wettability, a superhydrophilic taurine-electrografted cloth electrode was tested and found to result in a smaller medium frequency arc in the EIS data (∼0.4 Ω for all tested velocities, Fig. S7c), which validates our hypothesis of incomplete wetting for the studied electrodes. Nonetheless, even the performance of this superhydrophilic electrode (Fig. S8) cannot be completely captured with the model, which further reveals the complexity of simulating woven microstructures using simplified networks.
It is noted that an additional ohmic contribution is present in the system, as the experimental response can be captured using a correction procedure of k m and a resistance factor (Fig. 5d). This can be explained in the light of two factors: the network extraction and the distributed ohmic resistance within the system. As shown in the network extraction section, the inhomogeneous weave pattern of the cloth electrode cannot be precisely captured by pores and throats showing the limitation of the SNOW algorithm for woven materials. Furthermore, distributed ohmic resistance was proven by Pezeshki et al. to play a significant role in thick electrode structures. 18 This ohmic contribution is a result of the parallel combination of ionic and electronic resistances that occur in the liquid electrolyte and carbon fibers, respectively. This resistance is not captured by the high-frequency intercept in the EIS data, but by a distributed feature in the high to medium frequency range.
Despite the limitations of the SNOW algorithm to match the structural properties of the cloth electrode with complete accuracy, the modeling framework provides a versatile platform to analyze the impact of the computed network's structural properties on the overall cell performance from first principles. Further research should be conducted on using various geometrical shapes (e.g., Voronoi and Delaunay tessellations 66 ) to better represent the pores and throats in complex electrode geometries to increase the robustness of PNMs, as current network extraction methods are not optimal for multiscale images. 67 In addition, going forward, it will be important to ensure fully wetted electrodes in the experimental setup by performing e.g., thermal treatment, chemical etching, and chemical doping, to enhance the performance of the cells. Thus, the modeling of fully wetted systems is relevant, as incomplete wetting is undesired and should be overcome for practical applications. Nonetheless, incorporating the physics of multiphase flows to track gas pockets within the partially wetted porous structure is of interest and could be subjected to further research to increase the versatility and robustness of PNMs.
Porous electrode performance.-Pore network modeling can be used to understand the impact of the electrode microstructure on the performance of RFBs for distinct electrolytes, as local distributions of relevant properties (e.g., pressure, concentration, current) through the porous electrode can be obtained. This is an important advantage over experimental methods and macroscopic continuum models, as it enables spatially-resolved analysis of limiting phenomena; and over pore-scale direct numerical simulations, as it reduces computational requirements and affords simulation of larger electrode volumes. This provides valuable insight into the effects of microstructural properties on the electrode performance, useful for electrode optimization studies. Here, we compare the performance of both electrodes, assuming fully wetted systems, and demonstrate that the differences in macroscopic performance indicators directly correlate to the unique microstructural properties of both electrodes. To assess the performance on the cell level, the pressure drop, current density fields, concentration depletion effects, and overpotential contributions within both electrodes are discussed in the following subsections.
Pressure drop.-For efficient system operation, the pumping losses should be minimized compared to the electrical power output. Increased electrolyte flow rate enhances mass transport, improving the electrochemical performance, but simultaneously increasing pumping power requirements as the pumping power loss is defined as where P pump is the pumping power loss, ΔP the pressure drop, and η p the pumping efficiency. To keep the pressure drop minimal, the permeability of the electrode should be high, as the pressure drop and permeability are related by Darcy's law where L is the length of the medium and κ the permeability. The permeabilities of the extracted cloth were observed to be 14.4 (in-plane) and 17 (through-plane) times greater than those of the paper electrode, which resulted in a lower pressure drop (Figs. 6a-6b). The superior permeability of the cloth electrode is attributed to the larger pore segments (high porosity transport pathways) within the bimodal PSD (Fig. 3d).
Current density profiles.-The current distribution over the normalized thickness and length of both electrodes is shown in Figs. 6c-6f. Both electrodes show a distinct gradient in their current generation as a function of the electrode thickness, where most of the current is obtained in the region closest to the membrane. This distribution concurs well with the design principle behind zero-gap electrode cells, as the aggregated effect of ohmic losses is the smallest closest to the membrane. 68 This result highlights that the design of through-plane transport pathways of porous electrodes is not only important to reduce the pressure drop, but also to obtain efficient charge transfer within the porous electrode. Recently, it has been demonstrated that utilizing phase-separated induced porosity gradients in porous electrodes offers a feasible method to trade-off permeability and reactive surface area over the thickness of the electrode. 69 When operating at limiting-current conditions, the current distribution no longer displays a gradient as a function of the electrode thickness but shows a more uniform distribution throughout the electrode (Fig. S9) as all surface area is utilized for the redox reactions. Figure S9 thus shows uniform distributions at 1.5 cm s −1 for all systems except the cloth electrode with the nonaqueous electrolyte, which concurs well with Fig. 4 (no limitingcurrent conditions reached). Moreover, the aqueous electrolyte systems (Figs. 6c-6f) reveal a flatter drop-off in the current profile as a function of the electrode thickness, which can be explained by the higher ionic conductivity of the aqueous electrolyte (7.65 S m −1 ) compared to the non-aqueous electrolyte (1.99 S m −1 ). We anticipate that future work could leverage this systematic analysis to engineer the porous electrode thickness and local distributions of electrochemical surface area for specific chemistries and operating conditions. Reactant depletion.-Reactant depletion within the porous structure can be classified by longitudinal depletion and local depletion. The concentration profiles over the electrode length reveal that for both electrodes, the highest degree of reactant starvation occurs in the longitudinal dimension in the flow-through system (Fig. 7, at 1.5 cm s −1 for the non-aqueous electrolyte), as the inflow of fresh reactant is outweighed by the electrochemical consumption of reactant. At the end of the 1.7 cm electrode length (Figs. 7c, 7f), a slightly higher concentration is found for the cloth electrode in comparison with the paper electrode, which is a result of the thickness of the cloth electrode permitting a larger inlet flow rate at a constant inlet velocity.
The concentration distribution maps in Fig. 7 reveal that, while the bimodality of the cloth electrode was found to yield superior permeability, it results in a higher degree of local reactant depletion. Local clusters of smaller pores, with limited convective flow, show a higher degree of reactant depletion than the surrounding larger pore segments. Moreover, the unimodal PSD of the paper electrode results in a more uniform concentration profile.
Overpotential contributions.-To compare the electrochemical output of both electrodes in isolation, the performance must be corrected for the membrane resistance (see the calculations in section S5 of the Supporting Information) as different membranes were used for each electrolyte. The results are shown in Figs. 8a-8b, from which several interesting trends can be observed. First, the performance is similar within both microstructures for the nonaqueous electrolyte at high to moderate velocities; however, the cloth electrode outperforms the paper electrode at lower velocities due to its ability to reach higher limiting current densities. Second, the paper electrode outperforms the cloth electrode for the aqueous electrolyte at high to moderate velocities; however, the trend is reversed for lower velocities. Finally, the aqueous electrolyte outperforms the non-aqueous electrolyte, especially at high velocities. These trends can be further analyzed in the light of the individual overpotential contributions (i.e., activation, ohmic, and concentration) extracted from the symmetric flow cell simulations by the calculations presented in the Supporting Information (Section S5).
Two main trends can be observed in Figs. 8c-8f. First, due to the facile kinetics of both systems, the activation overpotential is the smallest contributor to the total overpotential (Figs. 8c-8d). The activation overpotential is inversely dependent on the exchange current density and thus the activation overpotential is higher for the aqueous electrolyte in both electrodes. Furthermore, as a consequence of the concentration dependency of the exchange current density, the activation overpotential increases with decreasing electrolyte velocities for all four systems.
Second, the concentration overpotential becomes dominant at low inlet velocities (Figs. 8e-8f). When comparing both electrodes, it is observed that the cloth electrode does not reach mass transfer limitations at the same current densities for both electrolytes. As explained in the reactant depletion section, care should be taken when attributing this effect to an improved mass transfer induced by the microstructure of the cloth electrode, as the obtained response is mostly a result of comparing the electrodes at the same inlet velocity. However, the pressure drop in the cloth electrode is still lower at similar velocities, indicating that the use of these electrodes is favored when operating in a high conversion regime. When comparing the performance at a velocity of 20 cm s −1 and −1 V, only 1-10% of the total overpotential is caused by depletion effects (Table SVI). Notably, the concentration overpotential at 20 cm s −1 is slightly higher in the cloth than in the paper electrode, which can be explained by reactant depletion effects in the highly electrochemically active regions consisting of smaller pores, as observed in Fig. 7. Furthermore, the concentration overpotential is greater for the aqueous electrolyte at a fixed potential, as the higher conductivity of the aqueous electrolyte enables operation at higher current densities, thus increasing the contribution of reactant depletion effects on the overpotential. In addition, the larger diffusivity of active species in the non-aqueous electrolyte, partly explained by the lower viscosity of acetonitrile, can help explain the reduced concentration overpotential for the non-aqueous electrolyte resulting in a faster rate of mass transfer of active species to the electrode surface.
In addition, for all chosen electrolyte-electrode pairs, the ohmic overpotential is dominant at moderate to high electrolyte velocities ( Fig. S10 and Table SVI), which can be explained by the facile kinetics, large membrane resistances, and the low effective conductivities of the chosen electrolytes (e.g., dilute conditions). The lower ohmic overpotential of the aqueous system can directly be explained by its higher effective conductivity.
The microstructural analysis showed that the bimodality of the porous cloth electrode had a positive effect on the overall permeation of fluids through porous electrodes, but the bimodal distribution should be engineered carefully to mitigate the effect of reactant depletion in highly electrochemically active regions. Furthermore, the balance between electrode thickness and surface porosity is critical as most of the current is generated near the membrane side of the porous electrode (when not operating at limiting-current conditions). This shows that the power output of the cell depends strongly on the electrode thickness as thicker electrodes will result in increased ohmic losses, but only in an incremental increase in the electrochemical output. To aid in the design of the porous electrode structure, bottom-up electrode optimization studies can be leveraged to balance these contradictory demands. Operation envelopes.-Besides the power of pore network modeling to investigate the impact of existing porous electrode microstructures, it can be leveraged to perform a broad screening of operating conditions and electrolyte chemistries for the identification of efficient operation envelopes. Here, we investigate the effect of the electrolyte velocity and redox chemistry on the cell performance in terms of electrical power output and pumping power losses (Fig. 9).
For efficient flow battery operation, two contradictory requirements must be met: (1) the electrical power output should be maximized, whilst (2) the pumping losses should be minimized. Therefore, we study the correlation between the pumping requirements and achievable current densities using the iR Ω -corrected electrochemical power. In this study, the iR Ω -corrected electrochemical power represents the useful electrical power driven by electrode performance and was calculated by where P el is the iR Ω -corrected electrochemical power, P theory is the theoretically achievable electrochemical power calculated by multiplying the open-circuit voltage (E OC ) of an iron-chromium battery of 1.18 V, chosen arbitrarily, by the current (I), and P losses are the power losses calculated by multiplying the iR Ω -corrected cell overpotential ( Ω E iR ) by the generated current. The pumping power was calculated using Eq. 31 and assuming a pumping power efficiency of 90%. 70 In Fig. 9, we plot the electrical power output vs required pumping power to identify operation envelopes that warrant efficient operation. In this figure, a velocity sweep was performed (at 5, 1, 0.5, 0.2, 0.05, 0.01, 0.005, 0.001, and 0.0005 cm s −1 ) for the aqueous electrolyte and two hypothetical electrolytes representing relevant systems (e.g., all-vanadium) with sluggish kinetics (exchange current densities of 1 × 10 −2 A m −2 and 1 × 10 −4 A m −2 , using the electrolyte properties of the aqueous electrolyte), for both electrode structures. The figure contains two demarcations indicated by the dashed lines illustrating a trade-off between the electrochemical power output and pumping power requirements 7 useful in the optimization of next-generation porous electrodes, electrolytes, and system parameters. We find that there is a plateau where higher electrolyte velocities only result in increased pumping requirements without achieving a significantly larger power output. Furthermore, we find that this plateau is reached at lower velocities for kineticallysluggish electrolytes and the higher permeability of the cloth enables a wider operational range. As a more extensive analysis, in section S6 in the Supplementary Information, the effects of various cell and electrolyte properties on the performance of the cell were further analyzed. We correlate the effect of the exchange current density and the inlet velocity to the overpotential distribution, outlet concentration, limiting current density, and pressure drop in both electrodes to show the potential of this modeling framework for systematic parameter screenings. Altogether, we have developed a pore network modeling framework and critically assessed its versatility and robustness. Based on our findings, and recommendations for further research on topological representations of woven electrodes and partial wetting, we believe this framework can be used as a predictive tool to engineer electrode microstructural properties, perform parameter sweeps, and investigate the effect of individual parameters on the system performance.

Conclusions
Engineering porous electrode microstructures in redox flow batteries is a promising approach to enhance stack performance and reduce pumping power requirements. To improve the fundamental understanding of the relationship between electrode microstructure and flow battery performance, we developed an electrolyte agnostic, microstructure-informed, and computationally inexpensive pore network modeling platform that solves for electrochemical kinetics, fluid flow, mass transport, and charge conservation. Compared to fully detailed pore-scale models (e.g., lattice Boltzmann simulations), the use of pore network modeling relies on geometrical simplifications which impact accuracy but enable simulation of larger sample sizes. Here, we demonstrated that the model can accurately describe electrochemical performance under a broad range of operating conditions and redox chemistries. Importantly, our systematic and rigorous analysis reveals some scenarios where the current modeling framework lacks accuracy (e.g., aqueous electrolytes with woven electrodes) which motivates further research on three-dimensional modeling of multi-scale microstructures and multi-phase flows.
The model was validated with experimental measurements with single electrolyte flow cells for the non-aqueous TEMPO · /TEMPO + electrolyte for two different electrode microstructures (Freudenberg H23 paper and the woven ELAT Cloth). We find that the experimental performance is well captured by the numerical model for this non-aqueous electrolyte without the use of fitting parameters. To further investigate the versatility and robustness of the pore network model, the performance of an aqueous Fe 2+ /Fe 3+ electrolyte was investigated within both electrodes. For this electrolyte, the simulations only partially describe the experimental performance. We hypothesize that incomplete wetting of the electrode structure with the aqueous electrolyte results in underutilized electrode surface area in the experimental measurements. Thus, we introduced an adaptation of the core physical descriptor of the mass transfer coefficient, resulting in a good fit for the paper electrode. Furthermore, we find that translating the woven electrode microstructure to a network connected by spheres and cylinders with pores spanning several length scales does not result in an accurate description of the electrochemical response. Thus, an additional resistance correction was necessary to obtain a better fit for the woven electrode with the aqueous electrolyte. Multiscale structures like carbon cloths motivate further research into alternative pore network geometries in addition to the investigation of incomplete wetting by experimental and modeling methods. Regardless, the model provides a robust platform to analyze the impact of the structural properties on the overall cell performance from first principles.
Moreover, the ability to simulate the pore-level performance of the electrode microstructure was discussed, showing the strength of using pore network models compared to experimental methods and continuum models, as spatially-resolved analysis of limiting phenomena can be performed, and direct numerical simulations, as larger electrode volumes can be simulated at a low computational cost. Here, it was observed that the presence of a bimodal structure of the woven electrode, with large pores within the threaded bundles, resulted in significantly reduced pumping power requirements in comparison to the unimodal pore size distribution of the paper electrode but also incited a higher degree of local areas with depleted reactant. The presence of local depletion zones explains the slightly larger concentration overpotential in the cloth compared to the paper electrode at moderate to high electrolyte velocities. Due to the low conductivity of both chosen model electrolytes, the ohmic overpotential was found to be the largest overpotential at velocities of 20 cm s −1 ; however, at velocities of 1.5 cm s −1 , longitudinal reactant depletion effects dominate the electrochemical losses at higher current density operation. Moreover, through-thickness current profiles show that most of the current is generated at the membrane-electrode interface at high electrolyte velocities, which motivates re-engineering electrode thicknesses for flow battery applications. We anticipate that future work could leverage this modeling framework to predict optimal electrode thickness and microstructure as a function of flow battery chemistry and operating conditions. Finally, we performed a parametric sweep to study the correlation between useful electrical power and pumping power requirements, hence identifying operation envelopes that warrant efficient operation. This knowledge can inform technology practitioners in selecting their porous electrodes for specific chemistries and to identify efficient operating conditions. In summary, the developed chemistry and microstructure agnostic pore network model provides a computational inexpensive framework to investigate porous electrode microstructures and allows for quick parameter screening over the electrode length. The model is well suited for system performance evaluations (electrolyte type, operating conditions), electrode microstructure comparisons, and analysis of flow field geometries in tandem with electrodes. Moreover, semi-quantitative modeling frameworks can enable the optimization of porous electrode structures for multiple system designs (e.g., redox chemistries, reactor architecture).
Going forward, the presented framework should be expanded by including complex flow field geometries, the modeling of the membrane 32 to simulate species crossover, and more detailed physics (e.g., migration, 30 ion-ion interactions, partial wetting) to better capture the electrochemical behavior of emerging nonsymmetric redox flow battery systems. Coupled with machine learning approaches, the proposed framework can be used as a predictive framework with the potential to accelerate the engineering of porous electrodes for specific electrolyte chemistries from the bottom-up.

Credit Statement
M.v.d.H. and R.v.G. contributed equally to this work. M.v.d.H. contributed to the conceptualization, methodology, experimental validation, network validation, operation envelopes study, formal analysis, investigation, data curation, writing-original draft, writingreview and editing, and visualization, R.v.G. developed the numerical model and contributed to the conceptualization, methodology, network validation, formal analysis, investigation, data curation, writing-original draft, writing-review and editing, and visualization. M.A.S. contributed to model review and debugging, and data curation. J.G. contributed to conceptualization, writing-review and editing, and supervision. Finally, A.F.C. contributed to the conceptualization, methodology, investigation, funding, resources, writing-original draft, writing-review and editing, project administration, and supervision.