Understanding the role of the porous electrode microstructure in redox flow battery performance using an experimentally validated 3D pore-scale lattice Boltzmann model

� A 3D pore-scale lattice Boltzmann model is developed to simulate flow and reactions. � Simulations are performed on 3D electrode microstructures from X-ray tomography. � The simulated performance is in excellent agreement with experimental data. � The model can be used to understand and design porous electrode structure. � Carbon cloth shows the desirable porous structure for high performance.


Introduction
Redox flow batteries (RFBs) are promising devices for large-scale energy storage due to long cycle life, high reliability, independent tunable power and energy, low cost and simplified manufacturing [1][2][3]. These batteries are used in applications such as load levelling and supplying reliable power and are well suited for storing electricity produced by intermittent renewable energy sources, such as wind or solar for long durations [3]. A typical RFB consists of two electrolyte reservoirs from which the electrolytes are circulated by pumping through two carbon-based porous electrodes. Electrochemical reduction and oxidation of redox couples occur at the negative electrode and the positive electrode, respectively, during charge, with the reverse processes occurring during discharge [4]. The electrodes in RFBs are responsible for providing active sites for redox reactions and facilitating the distribution of chemical species. Therefore, the performance of the RFB is dependent on the properties of the electrodes, in particular, their microstructure. State-of-the art, commercially available RFB electrodes consist of randomly distributed or aligned carbon fibers that form porous carbon mats with varying microstructural properties [4][5][6][7], and include non-woven structures (e.g., carbon papers and felts) and woven carbon cloths [7]. Although these materials are inert and durable, they also have some drawbacks including low surface area, poor wettability and high pressure drop [4,8]. To improve the battery performance, previous work has been dedicated to improving surface area and tuning surface chemistry [9][10][11]. For example, the incorporation of carbon nanotubes, carbon nanofibers or graphene oxide into the carbon electrodes have been used and proven to be effective in enhancing electrode performance [12][13][14]. Additionally, a number of surface modification methods, including acid [15,16] and thermal treatment [17,18], electrochemical activation [19,20] and other special surface treatments (e. g., plasma) [21,22], have been reported to enhance the electrochemical performance of electrodes. In these investigations, increased electrocatalytic activity is not only attributed to oxygen functional groups introduced onto the electrode surface, but also to the improved hydrophilicity [23].
Most approaches to date have targeted modification of electrode surface properties to improve the electrochemical kinetics. However, there has been less focus on understanding the role of the electrode microstructure on the performance of RFB. Recently, it was reported that introducing structural defects in electrodes through laser perforation for vanadium redox flow batteries (VRFBs) results in an up to 30% increase in power density [24,25]. Experimental studies on optimization of electrode microstructure for RFBs are very challenging as it is difficult to probe individual fibers or pores and fabricate different microstructures. The widespread availability of computational resources has made the use of advanced numerical modeling a valuable tool in assisting the design of electrodes to improve battery performance [26].
Numerical modeling can provide fundamental understanding of the coupled fluid flow, ion transport, electron conduction and electrochemical reactions that occur within RFB electrodes [26]. The numerical models that have been developed for the RFB can be classified into macroscopic continuum models and pore-scale models. In macroscopic continuum models [27][28][29][30][31], volumetric averaging conservation equations, such as mass, momentum, charge and energy, are solved in each elementary volume and the porous electrode is assumed as a homogeneous continuum with isotropic or anisotropic transport properties. Some key parameters, such as permeability of the porous electrode, ionic diffusivity through the membrane, the electrolyte viscosity and ionic interaction coefficients in the electrolyte affect the performance of RFBs, but are often determined using empirical relations which may not capture the real behavior due to the complex three dimensional structures of porous electrodes [27]. Most of the continuum models have been used to simulate cell performance. For example, a 2D continuum model for the VRFB was developed by You et al. [27] to investigate the effect of applied current density, electrode porosity and local mass transfer coefficient on the performance. Ma et al. [32] later extended that model to 3D using the same mathematical equations. The detailed spatial characteristics of concentration, velocity, overpotential and current density were obtainable with this model. Attempts have also been made to investigate the effect of the electrode properties on the cell performance. For example, a recent continuum model developed for a hydrogen-bromine RFB considered the effect of the porous electrode by giving different porosity and surface area, as well as the overall cell architecture [33]. Although useful for predicting the cell performance, the continuum models at the cell scale cannot be used to evaluate the precise effect of the electrode microstructure on the RFB performance. Alternatively, the pore-scale models represent an efficient way to obtain in-depth understanding of transport phenomena in porous media [34][35][36]. Qiu et al. [37,38] adopted the lattice Boltzmann method (LBM) to model electrolyte flow through a porous electrode, the structure of which was obtained using X-ray computed tomography (CT), while the finite volume method (FVM) was used to solve the coupled species and charge transport and predict the performance of the VRFB. Most recently, liquid-gas two-phase flow without reaction on the fibrous electrodes of VRFB was simulated using the LBM [39] to study the effects of porosity, fiber diameter, wettability of the fiber surface and saturation of the porous electrode on gas bubble migration and reduction of reactive surface area due to coverage of solid surface by bubbles. Thus, the 3D pore-scale LBM model provides a powerful tool to work with the real electrode microstructures and to investigate their effects on cell performance. Nonetheless, these models require experimental validation to verify the robustness of the conclusions drawn. However, direct comparisons between modeled performance and experimental measurement based on the same electrode materials applied in simulations and experiments have not been demonstrated in the literature to date.
In the current work, we investigate the influence of electrode microstructure conducted by combining computed tomography, electrochemical experiments, and 3D pore-scale modeling. Specifically, three different electrodes (Freudenberg paper, SGL paper, and carbon cloth) with different microstructures obtained by 3D X-ray CT were used for the LBM simulations of a non-aqueous redox flow cell. The experimental details are discussed in Section 2. The overpotentials for the three different electrodes were then calculated by solving the species/ charge/fluid transport coupled with electrochemical reaction on the surface of the carbon fibers, and the predicted results were validated against experiments (Section 3). Furthermore, the 3D distribution of current density and overpotential are presented and the effects of different microstructures and operation conditions are illustrated and compared with experimental data (Section 4).
In a redox flow cell operating in single electrolyte mode [40], the electrolyte is circulated through the high potential side, where the active species are oxidized. The outlet of the high potential side is connected to the inlet of the low potential side, where the active species are reduced before returning to the reservoir. An advantage of this setup is a constant state of charge, which makes this a platform well-suited for electrode comparison assuming other cell components are held constant. The details of the experimental setup and the assembly of the cells can be found in our paper [40].
Three commercially sourced carbon electrodes (details are given in Section 2.2) were employed in the single electrolyte setup to test the electrochemical performance. The polarization measurements were performed using an Arbin battery tester (FBTS-8) by potentiostatic holds of 30 s in 25 mV steps between 0 and 0.6 V, recording one data point per second. The mean of the current and the voltage were calculated for the final 50% of data points to ensure steady state conditions. The ohmic resistance of cells for iR-corrected polarization curves was determined from the high frequency intercept of Nyquist plots obtained from electrochemical impedance spectroscopy [40].
Pressure drop through different electrodes was measured under relevant flow cell operation conditions in a custom half flow cell setup in which the separator is replaced with a solid copper plate (3 mm thickness, McMaster-Carr) to prevent leakage, described in more detail in Ref. [40]. Measurement of the pressure drop was done using digital pressure gauges (SSI Technologies). Pressure loss was obtained as a pressure differential of outlet and inlet pressure, and the pressure loss through an empty cell was subtracted to isolate the porous electrode contribution. MeCN without active species was used as the fluid for pressure drop experiments rather than redox electrolytes since both fluids are Newtonian with similar surface tension and viscosity, and therefore have similar wetting properties.

Tomographic imaging and data processing
The microstructure of three commercially sourced carbon electrodes (Freudenberg H23 paper from Fuel Cell Store, SGL 399 AA paper from Ion Power Inc. and Nuvant ELAT carbon cloth from Fuel Cell Etc) were imaged using a laboratory X-ray CT system (Phoenix Nanotom 180, GE, USA). Table 1 shows the tomographic imaging experiment parameters employed for each electrode sample.
At the respective X-ray tube voltage and current settings for each sample, several projection images were captured over a 360 � sample rotation. Tomographic reconstruction of the acquired projection images was performed with Datos-x reconstruction software (GE, USA), which uses a cone-beam filtered back-projection algorithm based on the Feldkamp-Davis-Kress (FDK) algorithm [43]. Image pre-processing and segmentation of the resulting reconstructed volume was performed using commercial image processing software (Avizo v9.5, Thermo Fisher Scientific). A 3D median smoothing filter was applied to the reconstructed volume to reduce the image noise. Threshold-based segmentation of the volume was then performed to isolate the carbon fibers from the pore/electrolyte phase.

Fluid transport
A three-dimensional single-relaxation-time (SRT) LBM is adopted to simulate the liquid electrolyte through the connected pore space in the porous electrode. In the LBM model, the motion of a fluid is described by a set of discrete single-particle density distribution functions. The particle distribution function with a single relaxation time can be written as: where f α is the particle distribution along the α th direction, f eq α is the equilibrium distribution, δt is the time step, e α is the particle velocity in the α th direction and τ is the single relaxation time. The equilibrium distribution function f eq α ðx; tÞ [44] in Eq. (1) can be calculated as: where w α is the weighting factor, c s is the lattice sound speed and v is the velocity of fluid. In this study, the D3Q19 model is adopted for the fluid transport simulations. The discrete velocities and weighting factors for D3Q19 are given by Ref. [44]: The corresponding macroscopic density and velocity are calculated by: e α f α ðx; tÞ (6) where N is the number of discrete particle velocities.

Simulation of species transport
In the present study, only the transport of the active species j 2 fTEMPO; TEMPO þ g is considered due to the nature of the single electrolyte experiment. The concentration of BF 4 is obtained via the electroneutrality condition, P z j C j ¼ 0. The transport of species j within the electrolyte is governed by the convection-diffusion equation [45]: where C j is the concentration of species j, D j is the diffusivity of species j, z j is the charge of species j, R is the universal gas constant, T is the absolute temperature (K), and ϕ is the electrical potential. r⋅ðC j vÞ is the convection term and r⋅ðD j rC j Þ is the diffusion term. r⋅ where c j;α is the concentration distribution function, e α is the particle velocity in the α th direction, and τ Dj is the dimensionless relaxation time used to determine diffusion coefficient. For porous media with high porosity (>80%) such as the carbon paper electrodes used in this work, the D3Q7 lattice model without lattice velocity i ¼ 7-18 is sufficient to accurately simulate species transport in electrolyte [35,39]. The equilibrium function is specified as: where C j is the concentration of the species j which is obtained by: and A α in Eq. (9) is given by Ref. [47]: where A 0 is equal to 0.25 and D j ¼ ð1 A 0 Þðτ Dj 0:5Þ= 3. The Neumann boundary condition is applied as the wall concentration boundary condition which can be described as D j j is the reaction rate per unit area for different species (a term explained in greater depth in 3.1.3) and n is the local surface normal pointing from the solid phase towards the electrolyte phase [45]. To implement this boundary condition with LBM, the distribution function entering from the wall c j;α ðx; t þδtÞ should be determined as [48]:

Electrochemical reaction
A single redox reaction takes place within the single electrolyte flow cell occurring at the solid-liquid interfaces: The molar flux for TEMPO and TEMPO þ at the positive electrode during charge will be [49]: where k is the reaction rate constant, F is the Faraday constant and α is the transfer coefficient. In the current LBM model, C S j is the concentration of species j on the surface of carbon fibers, while C j is the bulk concentration of species j in the pore near the carbon fibers, with j being TEMPO or TEMPO þ . The consideration of the concentration of species j on the surface of carbon fibers (C S j ) and the bulk concentration of species j in the porous space near carbon fibers (C j ), is to include the mass transfer loss by the species transport between the bulk solution in the pore and the carbon fiber surface, which has been confirmed to be important for RFBs [49]. The difference between the bulk concentration in the pore and the surface concentration of species ΔC j , is related to the reaction rate r '' j , by introducing the mass transfer coefficient k m [50]: k m , which characterizes the mass transfer loss, is dependent on the flow rate and velocity of the electrolyte in the pore space. In this study, the correlation between the mass transfer coefficient k m and the velocity of the electrolyte is obtained by fitting the simulated results with the experimental results for three different electrodes. An empirical equation to express the correlation between k m and the velocity is derived through this fitting process: The velocity v used in Eq. (16) is the local velocity of electrolyte near the carbon fiber, and is decided by the flow rate. The mass transfer coefficient is proportional to the electrolyte velocity near the carbon fiber surface and tends to become a constant value as the velocity of electrolyte continues to increase. The overpotential η at the electrodeelectrolyte interface for the single electron transfer is defined by the Nernst equation [49,50]: where ϕ s and ϕ e are local potentials for electrode and electrolyte, respectively, and E 0 is the standard equilibrium potential electrode.

Simulation of charge transport
Inside the solid carbon fibers, the governing equation for the potential field can be written as: where κ s is the electrical conductivity of the solid fibers.
In the electrolyte, the potential field can be obtained as: where the effective electrical conductivity of electrolyte κ e is defined as [37,38]: The LBM is used to solve Eq. (19) and Eq. (20). The LBM equation is written as follows: where W is equal to zero for Eq. (19). The Neumann boundary condition is applied at the electrode-electrolyte active interface: where J is the current density that can be calculated by Faraday's law: To implement this boundary condition with LBM, the distribution function entering from wall g α ðx; t þδtÞ should be determined as [48]:

Calculation of the cell concentration overpotential and activation overpotential
The local concentration overpotential within the electrode microstructure is calculated from the Nernst equation as follows [50]: If the difference between the concentration in the pore and the concentration at the fiber surface is negligible (C i ¼ C S i ), then the concentration overpotential term in Eq. (26) will also become negligible. The average concentration overpotential over the 3D electrode microstructure consisting of carbon fibers from the Nernst equation is obtained by integrating local concentration overpotential η conc Nernst , over the surface of the solid fibers η conc;Cell where S CC is the surface area of the current collector. Under conditions of negligible mass transfer resistance C i ¼ C S i , the voltage losses at the electrodes are purely due to charge transfer resistance. In that case, the concentration overpotential from the Butler-Volmer equation will disappear. Eq. (14) is reduced to Ref. [50]: The η act represents the voltage losses to overcome the activation barrier, and is, therefore, referred to as the activation overpotential. So, the activation overpotential of the cell is calculated by integrating η act over the solid fiber surface as follows: The concentration overpotential component arising from the Butler-Volmer equation can be calculated by: The total concentration overpotential of the cell is the sum of both the Nernst (thermodynamic) and the Butler-Volmer equation (kinetic) contributions: η conc;Cell ¼ η conc;Cell BV þ η conc;Cell Nernst (31) Note that we introduce the mass transfer coefficient k m as in Eq. (16), to obtain good fit between the experimentally measured and simulated overpotentials. Reducing the value of k m would result in overpredicted concentration overpotentials. On the contrary, increasing k m results in underpredicted concentration overpotentials. Fig. 1 shows the 3D porous electrodes reconstructed from X-ray CT. The inherently different microstructure of the three electrodes applied in the RFB cell is obvious from Fig. 1. The Freudenberg paper has densely and randomly packed fibers, yielding narrowly distributed pore sizes with a dominant peak around 20 μm, with all the pores below 60 μm. The SGL paper consists of randomly distributed fibers intersecting with each other, giving a wide pore size distribution in the range of 0-160 μm, with two dominant peaks around 10 μm and 100 μm, and a big fraction of pores (more than 50%) between 60 μm and 160 μm.

Computational details
Carbon cloth, unlike the Freudenberg paper and the SGL paper, has a woven and ordered microstructure from fiber bundles, giving a dominant pore size peak around 10 μm (arising from inside the fiber bundles) and a small fraction (less than 20%) of bigger pores spanning from 60 to 160 μm (arising from between the bundles).
Note that fresh samples without compression were used in X-ray imaging. For experimental performance testing, the samples were compressed at the point of battery assembly, with a compression ratio of 10-20%. The compression is applied along the electrode thickness direction to increase the fiber contact and the electronic conductivity. We expect the porosity of the compressed samples in the experiment is not significantly different from that of the fresh samples. According to Jervis et al. [51], when the compression ratio is less than 25%, the compressed electrode microstructure remains almost the same to the fresh samples, with less than 5% change in porosity and negligible change in surface area. Future work should study the influence of compression on mass transport overpotential. Fig. 2 depicts a single electrolyte flow cell set up for the simulation work, which shows the liquid electrolyte flow from the negative electrode (anode) to the positive electrode (cathode), whilst the two electrodes are separated by the membrane. For each type of the electrode material, the same microstructure, extracted from the tomography data, is used for the anode and the cathode. This is done to mimic the experimental setup.
In this paper, it is assumed that the electrode is completely saturated with electrolyte and all the fiber surfaces within the electrode participate in the reaction. This assumption is reasonable considering the low contact angle of acetonitrile on the carbon surfaces, which is observed in our experiments [40]. It is noteworthy that the LBM model can also simulate the gas-liquid two phase flow within the 3D porous electrode. In fact, our previous paper has demonstrated that in some cases, gas bubbles could be trapped within the porous structure, leading to reduced available surface area for reactions [52]. The detailed boundary conditions for species and charge transport are shown in Table 2.
The membrane is not included in current simulations as the crossover effects do not degrade the performance of the single electrolyte flow cell [42]. Migration of species across the membrane is neglected and a zero-gradient flux boundary condition is imposed for all species. The half bounce back scheme in LBM is applied on the wall boundary in this study. Constant velocity boundary conditions are imposed on the inlet to control the flow rate of electrolyte. Table 3 illustrates the physical parameters used in the simulation. Table 4 shows the sizes (i.e. the computational domains) of different electrodes used in LBM simulations.
The electrode structures used in the LBM simulations share the same thickness (T) with the actual electrodes used in the experiment, while smaller length (L) and width (W) are adopted to reduce the computational cost. Because of the reduced cell active area in simulations, the flow rate values needed to be scaled down to ensure they had the same ratio between reactants supply rate and consumption rate when the reactions in the RFB achieved an equilibrium state. Here, a dimensionless number R is introduced to describe the ratio between the reactant supply rate and consumption rate, which is defined as: where J is current density, Q is the flow rate which is calculated as Q ¼ V average A inlet and S is the cell active area S ¼ WL. In the current work, the dimensionless number R is adopted to normalize the flow rate for different cases in experiments and simulations with the same current density. Recognizing that the simulated electrode structures (as given in Table 4) are smaller than a typical RFB electrode, we want to test if the simulated structures are the representative elementary volume. This is done by comparing the properties of the simulated structures (which is a subset of the actual electrode) and the actual electrodes (denoted as "Master"). As shown in Table 5, porosity and specific surface area are selected as points of comparison between the subset and the master structure. Very similar porosity and specific surface area values are obtained for both the subset and the master structures, for all the three different electrode materials. In addition, the relationship between pressure drop per unit length and flow rate from three different simulated structures are compared with the results from experiments, as shown in Fig. 3. Excellent agreement is observed between the simulated and the experimental pressure drop. As the pressure drop of fluid flow in porous media is sensitive to the geometrical properties of a porous structure [56,57], we conclude that the chosen computational domains are sufficient to represent the three different electrodes. As can be seen in Fig. 3, different electrode materials give different pressure drop behaviors. The Freudenberg paper shows a much higher pressure drop than the SGL paper and carbon cloth. This is because the Freudenberg paper has a densely packed structure and a narrow pore size distribution around 20 μm (as shown in Fig. 1), which prohibits the flow of the electrolyte through the structure. The presence of larger pores in both the SGL paper (with pore size up to 100 μm) and carbon cloth (with pore size up to 160 μm) is beneficial for the flow of the electrolyte, and thus the reduction of the pressure drop in the structure. The presence of pores in the range of 100-160 μm in carbon cloth gives a lower pressure drop compared to SGL paper, confirming again the benefit of having large pores, even when such large pores only take up a small fraction in the whole pore size distribution (as shown in Fig. 1).
To check the effects of the grid size on the LBM model results, simulations for the Freudenberg paper and carbon cloth are carried out for four different grid sizes respectively. The results of the iR-corrected overpotential (η Cell act þ η conc;Cell BV þ η conc;Cell Nernst ) for different grid sizes are presented in Table 6. The values of the iR-corrected overpotentials are nearly the same for the grid resolutions of 3 μm, 2 μm and 1.5 μm, with a maximum overpotential difference of 2%. However, the over-potential for the grid size with l LB ¼ 4 μm has 10% higher than the result from l LB ¼ 3 μm. Therefore, the grid size of 3 μm is chosen in this study for the simulation efficiency and accuracy.

Results and discussions
The numerical results are first compared with experimental results to validate the model. In the experimental study, it is difficult to distinguish the concentration overpotential and activation overpotential. Electrochemical impedance spectroscopy responses are complex and it is difficult to deconvolute overpotential contributions without a physicsbased model. In our LBM simulations, the value of concentration overpotential and activation overpotential can be calculated separately. The effects of the operating parameters and porous structure of electrodes on the concentration overpotential and activation overpotentials are investigated according to the results from LBM simulations. Fig. 4 shows the simulated and experimental values of the iRcorrected overpotential of the redox flow cell at different current densities and flow rates for three different electrodes (Freudenberg paper, SGL paper, and carbon cloth). The iR-corrected overpotential is the sum of activation overpotential and concentration overpotential. It can be seen from Fig. 4 that the simulated results are in good agreement with the experimental data over a wide range of current densities and flow rates for all three electrodes. This indicates that the 3D pore-scale LBM model is able to give reliable simulations of the performance of the RFB electrodes, and thus can be used to predict the RFB electrode performance. The iR-corrected overpotential increases as the current density increases and decreases as the electrolyte flow rate increases. At high flow rates (R > 13.8), the iR-corrected overpotential increases linearly with the increase of the current density. At low flow rate (R < 13.8), the iR-corrected overpotential increases exponentially with the increase of the current density. From Fig. 4, different iR-corrected values are observed for different electrodes. The effects of different electrodes on the electrochemical performance are further investigated and discussed in the following sections.

Investigation of concentration overpotentials
The 3D pore-scale LBM model is able to simulate concentration profiles within the electrode microstructure and derive the concentration overpotential. Fig. 5 (a) and Fig. 5 (b) show a 2D local distribution of concentration and the corresponding averaged profiles as a function of R and current density, for the Freudenberg paper electrode. The nondimensional number R, defined in Eq. (32), represents the electrolyte flow rate. In this study, the electrolyte [TEMPO, TEMPO þ ] flows into the pore space of the electrode at a concentration of 250 mol m 3 from the inlet as shown in Fig. 5 (a). The electrolyte from the positive side flows directly into the negative side in order to be consistent with the experimental setup shown in Fig. 2. The concentration near the inlet  approaches 250 mol m 3 as the Dirichlet concentration boundary condition is used. As shown in Fig. 5 (b), the concentration of TEMPO decreases along the length of the positive half-cell as the electrochemical reaction is converting TEMPO to TEMPO þ , while the amount of TEMPO increases on the negative side due to the reverse reaction. At low flow rate (i.e., R ¼ 4.13), significant polarization between the positive and the negative sides is observed. Compared to the cases of high flow rate (i. e., R ¼ 55.2), the low flow rate case (i.e., R ¼ 4.13) shows a much bigger difference in the concentration profiles between the positive and the negative electrodes. In Fig. 5 (b), the minimum average concentration of TEMPO on the positive side approaches 30 mol m 3 and the maximum average value on the negative side is close to 290 mol m 3 for the case with R ¼ 4.13. However, for the high values of R ¼ 55.2, the concentration profiles on the positive and the negative sides are almost symmetric, indicating low polarization. The minimum average concentration of TEMPO on the positive side approaches 220 mol m 3 and the maximum average value on the negative side is close to 280 mol m 3 . At the same flow rate value of R ¼ 55.2, higher current density (J ¼ 408.7 mA cm 2 ) gives slightly higher polarization of concentration, as shown in Fig. 5 (a). Higher current density means faster reaction rate (Eq. (24)), resulting in faster consumption of reactants.
We expect that such effect of R and current density on the concentration profiles will also be observed for the concentration overpotentials. This is demonstrated in Fig. 5 (c) and Fig. 5 (d). Fig. 5 (c) shows the 3D distribution of the concentration overpotential at the surface of Freudenberg paper electrode on the positive side for a small part (from x ¼ 45 μm to x ¼ 81 μm) of the simulated porous electrode with different values of R and current density. With the decrease of R, we see the increase of concentration overpotential on the positive electrode, due to the lack of electrolyte supply (low flow rate). Fig. 5 (d) illustrates the average concentration overpotential along the length of the cell (Y direction in simulation). The average concentration overpotential is calculated by integrating the local concentration overpotential on all of the carbon fiber surfaces. For the case of R ¼ 4.13, the average concentration overpotential on the positive side reaches its maximum value (�350 mV) near the outlet (Y � 300 μm), then drops to nearly 40 mV near the outlet on the negative side. In contrast, the concentration overpotentials are quite small with R ¼ 55.2. At the same value of R ¼ 55.2, the increase of concentration overpotential can be observed with the increase of current density due to the increased electrolyte consumption.
In order to compare the effect of different carbon electrodes on the performance, the concentration overpotential for the three different electrodes at the same flow rate and different current densities are compared in Fig. 6. At the same operation conditions, the SGL paper gives much higher concentration overpotential, compared to carbon cloth and the Freudenberg paper. The higher concentration overpotential in the SGL paper electrode could be attributed to the nonuniform and widespread pore size distribution of the SGL paper, as compared to the relatively uniform and narrow pore size distribution of the Freudenberg paper and carbon cloth, as shown in Fig. 1.
The pore size distribution impacts the electrolyte velocity field within the porous structure, resulting in significant differences in performance for each of the electrodes. Fig. 7 shows the velocity field and the TEMPO concentration distribution in these three different electrodes. In the Freudenberg paper electrode, most of the pores are concentrated around 20 μm in diameter, showing a narrow pore size distribution, so that it has relatively uniform flow velocity distribution ( Fig. 7 (a)) and small variance in concentration distribution inside the electrode ( Fig. 7 (b)). The carbon cloth electrode shows a dominant peak of pore size distribution around 10 μm, with a small amount of bigger pores in the range of 10-160 μm, which is different from the Freudenberg paper. Although it has a similar concentration overpotential compared with the Freudenberg paper electrode, the pressure drop in the carbon cloth electrode is much lower than in the Freudenberg paper (as shown in Fig. 3), due to the presence of larger pores. Compared to Freudenberg paper and carbon cloth, SGL paper does not show dominant peaks. Rather, it shows a non-uniform distribution over a wide pore size range of 0-160 μm, which leads to the non-uniform flow distribution within the electrode (excluding the region near the inlet and outlet). In that case, as shown in Fig. 7 (a), some regions of porous electrode are quite difficult to access by the electrolyte reactant, while most of electrolyte will select the path of the least resistance through the largest pores of the electrode. It can be seen from Fig. 7 (b) that higher concentration polarization occurs in some regions of electrode, which results in a higher concentration overpotential in the SGL paper electrode.

Investigation of activation overpotentials
In addition to the concentration overpotential, the activation overpotential also contributes to the total overpotential of the redox flow cell. Fig. 8 (a) shows the effect of flow rate on the average activation overpotential at the same external current density (J ¼ 408.7 mA cm 2 ) and different flow rates. From Fig. 8 (a), larger variations in both the average activation overpotential and local current density along the Ydirection are seen at the low flow rate of R ¼ 4.13 compared to the case  at the high flow rate of R ¼ 55.18. The maximum value of the local activation overpotential in the Freudenberg paper electrode approaches 6 mV when R ¼ 4.13, which is higher than the maximum local activation overpotential of 2 mV when R ¼ 55. 16. When current densities are the same (J ¼ 408.7 mA cm 2 ), we see larger variation in the concentration profile at low flow rate (R ¼ 4.13), compared to at high flow rate (R ¼ 55.16). Therefore, we can postulate that the large variation in activation overpotential at low flow rate, as observed from Fig. 8, is closely associated with the large variation in concentration distribution. This is consistent with Eq. (14) and Eq. (24), in which a large variation in concentration distribution will cause large activation overpotential when current density is kept constant. Fig. 8 (b) shows the average activation overpotential and current density along the Y-direction with the same flow rate and different external current densities for the Freudenberg paper electrode. The external current density is the current density at the current collector. The average activation overpotential at large external current density (J ¼ 408.7 mA cm 2 ) is higher than that at small external current density (J ¼ 205.9 mA cm 2 ) at the same flow rate (R ¼ 55.2). Fig. 8 (b) shows the average activation overpotential increases nearly by two times when the external current density increases from 205.9 to 408.7 mA cm 2 . This is because the magnitude of average activation overpotential is proportional to the value of the external current density. The local current density is obtained directly from the Butler-Volmer equation (Eq. (14)) and Faraday's law (Eq. (24)), in which the local current density is proportional to the local activation overpotential, so that the trend in the local current density is identical in shape to the trend observed for local activation overpotential on the surface of fibers.
In order to evaluate the effect of the porous structure on the activation overpotential, the activation overpotential from carbon cloth, Freudenberg paper and SGL paper electrodes are calculated from simulations at the same value of flow rate and external current density. As the difference in thickness will lead to the variation in the flow rate in an electrode with the same velocity, a same thickness (Z ¼ 156 μm) is set up for the three different electrodes in this particular simulation work. Fig. 8 (c) shows the overall activation overpotential values for the three different electrodes at the same flow rate. From Fig. 8 (c), we can see that all the three carbon electrodes give low activation overpotentials (below 10 mV) for the range of current densities investigated. These values are notably much lower (one order of magnitude lower) than concentration overpotentials reported in Fig. 6. This is attributed to the fast kinetics of TEMPO/TEMPO þ redox reaction on the carbon surface, as reflected by the high reaction rate constant k (obtained from the experimental measurement) in Table 3. For other redox reactions, higher activation overpotentials may be expected.
It is obvious that the SGL paper electrode has higher activation overpotential, compared to carbon cloth and Freudenberg paper electrodes. This could be attributed to the larger variation in concentration distribution inside the SGL paper, as seen in Fig. 7. As can be seen from Fig. 8 (c), there is a non-negligible difference in the activation  overpotential for Freudenberg paper and carbon cloth, with the Freudenberg paper giving the lowest activation overpotential. This is different from the concentration overpotential as shown in Fig. 6, where the Freudenberg paper and the carbon cloth display a similar performance. This indicates that other factors, apart from the concentration distribution, also affect the activation overpotential. One should also be aware that, the local distribution of activation overpotential and current density are dependent on the physical distribution and local connectivity of the carbon fibers. The intrinsic properties of porous electrode, such as the specific surface area, the connection of fibers, and the fiber surface properties, will also lead to the variation in activation overpotential in different electrodes. Future work should examine all the physical and chemical properties of different electrodes to fully understand the effect of the electrode on activation overpotentials.

Conclusions
A 3D pore-scale LBM model is developed to simulate electrolyte flow, species and charge transport as well as electrochemical reaction in a non-aqueous redox flow cell, serving as a platform to investigate the performance of different electrode microstructures (Freudenberg paper, SGL paper and carbon cloth). The detailed structures of the electrodes are obtained using X-ray computed tomography, and are input into the LBM model for simulating the electrochemical performance. A single electrolyte flow cell is employed to measure the iR-corrected overpotential (which is the sum of concentration overpotential and activation overpotential) for the different porous electrode structures as a function of electrolyte flow rate and current density. The simulated iRcorrected overpotentials are compared with the results of the single electrolyte experiment and showed excellent agreement over a wide range of flow rates and current densities. The values of concentration overpotential and activation overpotential obtained from LBM simulations are found to be dependent on the concentration distribution inside the electrode. Higher concentration and activation overpotentials are observed with large concentration polarization inside the electrode.
Furthermore, it was observed that electrode microstructure plays an important role on the concentration overpotential and the activation overpotential. The SGL paper electrode, which features a broad pore size distribution with multiple dominant pore size peaks (mainly around 10 and 100 μm) over a wide range of pore sizes (0-160 μm) and a big fraction of pores (more than 50%) in the range of 60-160 μm, results in a higher concentration overpotential. This is in contrast with the Freudenberg paper which features a single domain peak around 20 μm and no pores in the range of 60-160 μm, and carbon cloth which features a single dominant peak around 10 μm and a small fraction (less than 20%) of pores in the range of 60-160 μm. The non-uniform flow distribution inside the SGL paper electrode leads to large variations in the concentration distribution and therefore higher overpotential. The   Fig. 8. (a) Average activation overpotential, current density and concentration of TEMPO along Y-direction with same external density (J ¼ 408.7 mA cm 2 ) and different flow rates in Freudenberg paper electrode; (b) Average activation overpotential and current density along Y-direction at the same flow rate (R ¼ 55.2) and two different external current densities in Freudenberg paper electrode; (c) Activation overpotential of cell for three different electrodes at the same flow rate (R ¼ 4.13).
concentration distribution also affects the activation overpotential, and the SGL paper shows the highest activation overpotential among the three electrodes. The results also show that the activation overpotentials for the Freudenberg paper and carbon cloth do not follow the trend in their concentration overpotentials. This indicates that other factors may come into play and that further studies are required to fully understand the effect of the electrode microstructure on activation overpotential. The trend we observed with the impact of different porous structure of the electrodes on the concentration overpotentials should be applicable for other redox pairs. However, with kinetically-sluggish systems (i.e., all-vanadium RFB), the role of surface chemistry and electrochemically active surface area are expected to play a major role, since both influence exchange current density and, thus, activation overpotentials. In this case, our learnings (i.e., woven electrode provide favorable mass transfer) will still apply, but additional kinetic considerations are required (e.g., surface treatment to increase surface area, functionalization to increase kinetic activity, etc.). In other words, the woven electrode microstructure, providing optimal mass transport, can be used as a scaffold material in tandem with advanced surface modifications, to enable optimal kinetic performance.
Pressure drop through the three electrodes show that the Freudenberg paper (with a narrow pore size distribution around 20 μm and all the pores below 60 μm) has the highest pressure drop, whilst carbon cloth (with a dominant pore size peak around 10 μm and a small fraction of large pores in the range of 60-160 μm) has the lowest pressure drop.
The results confirm that the presence of large pores is beneficial for the electrolyte flow and the reduction of pressure drop. Considering that a high pressure drop is associated with a high operational cost, a porous structure with large pores may be desirable for maintaining low operational cost. Furthermore, a small fraction of large pores in the pore size distribution may be sufficient to have low pressure drop. From this study, we can conclude that, to have both good electrochemical performance and low pressure drop, a desirable porous structure should have a single dominant pore size peak around 10-20 μm and a small fraction of large pores. Overall, the 3D pore-scale LBM simulations predict cell performance with good accuracy, and it can be used to elucidate structure-performance relationships for RFB electrodes. This platform can be then leveraged to design optimized electrode architectures for high-performing redox flow batteries.