Exploring the Impact of Electrode Microstructure on Redox Flow Battery Performance Using a Multiphysics Pore Network Model

The redox flow battery is a promising energy storage technology for managing the inherent uncertainty of renewable energy sources. At present, however, they are too expensive and thus economically unattractive. Optimizing flow batteries is thus an active area of research, with the aim of reducing cost by maximizing performance. This work addresses microstructural electrode optimizations by providing a modeling framework based on pore-networks to study the multiphysics involved in a flow battery, with a specific focus on pore-scale structure and its impact on transport processes. The proposed pore network approach was extremely cheap in computation cost (compared to direct numerical simulation) and therefore was used for parametric sweeps to search for optimum electrode structures in a reasonable time. It was found that that increasing porosity generally helps performance by increasing the permeability and flow rate at a given pressure drop, despite reducing reactive surface area per unit volume. As a more nuanced structural study, it was found that aligning fibers in the direction of flow helps performance by increasing permeability but showed diminishing returns beyond slight alignment. The proposed model was demonstrated in the context of a hydrogen bromine flow battery but could be applied to any system of interest. © The Author(s) 2019. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0721910jes]

Solar and wind are the most prominent renewable energy sources and are increasingly being deployed at a surprising growth rate of 3% per year. 1 However, these sources are intermittent and it has been shown that consequently, integration of more than 20% renewable energies in the current aging grid endangers its stability. 2 Furthermore, renewable energy sources are active mostly during daytime whereas power consumption peaks during the evenings. Therefore, the generated power from renewable energy sources is often wasted due to lack of customer demand at the time of generation. To this end, energy storage technologies are actively being considered to alleviate the unpredictability and variability of renewable energy sources. [3][4][5] During the past decades, much of the focus of energy storage research has been in the areas of transportation and portable power where space/volume requirements are the limiting design factors, leading to advanced Li-ion batteries and hydrogen fuel cells. 6 However, these requirements are less critical in grid storage applications. Moreover, large-scale application of batteries in grid storage necessitates a reasonable number of charge/discharge cycles, ability to scale-up, and acceptable capital costs. Redox flow battery (RFB) technology has been demonstrated as a promising solution that meets many of these requisites. 7 The RFB can be viewed as a conventional rechargeable battery except the fuel is decoupled from the reaction sites and is pumped through the electrodes on demand. Similarly, a typical RFB cell consists of a cathode, an anode, and a membrane that connects the two electrodes and ensures electroneutrality while preventing cross-mixing.
Despite the advances in the RFB technology, they are not economically feasible yet and require further optimization for large-scale commercial deployment. Major losses in the RFB have been identified to be those related to mass transport limitations, rather than sluggish kinetics. Having said that, finding new redox pairs and catalysts is still an active area of research. [8][9][10][11][12][13] Unfortunately, since typical RFB electrodes are only a few hundred microns thick, the local concentration of species and the electric potential of the electrolyte cannot be reliably measured. To this end, modeling has proven to be useful in improving the electrode design and eventually to minimize the mass transfer losses. There have been numerous and extensive studies on modeling the RFB performance, ranging from control-oriented [14][15][16] to numerical models including mostly continuum [17][18][19][20][21][22][23][24] and only recently detailed pore-scale models. [25][26][27] In the modeling literature to date, the vast majority of papers focus on developing accurate continuum macroscopic models, which by their very nature are insensitive to the microstructure of the RFB electrodes. Though these models are helpful in predicting the effect of macroscopic features such as flow configuration, channel dimensions, and operating conditions on the RFB performance, they cannot be used to predict how the microscopic features such as the internal structure of the electrode affects the RFB performance.  for the first time introduced a multiphysics porescale model to study a vanadium RFB. 28 They used the lattice Boltzmann method (LBM) to obtain the velocity field in a digitally reconstructed electrode, and eventually used a finite volume-based solver to obtain the local concentrations of species and the electric potential of the electrolyte. In follow-up work, they used the model to study how electrolyte flow rate and electrode porosity affects the overall performance of the RFB. 29 Recently, Chen et al. (2017) published a similar study with a few differences. 27 First, they used LBM to solve both momentum and mass transport equations. Second, they did not include electrochemistry as a separate coupled physics, but rather assumed a constant fixed overpotential throughout their simulations. The major contribution of their work, however, is that for the first time they considered bubble formation in the electrode caused by oxygen evolution side reaction. This is important since it directly impacts macroscopic properties such as the effective diffusivity of species and can potentially be used for upscaling to continuum models. Finally, Fetyan and Banerjee et al. very recently published two on pore-scale characterization of carbon-carbon composite electrospun RFB electrodes in which they used pore-network modeling to compute effective transport properties of the electrode. 30,31 However, their study is focused on characterization of RFB electrodes rather than investigating the pore-scale effects of overall performance. A similar study has also been reported by Kok et al. but is also limited to computing effective mass transfer coefficients of flow battery electrodes using LBM simulations. 32 To the best of our knowledge, no other original work has been published on pore-scale simulation of the RFB, while continuum-based modeling and optimization of the RFBs are still being actively pursued. [33][34][35][36] This gap in the modeling literature is partly attributed to the lack of a general pore-scale modeling framework with feasible computing requirements. Both the works of  and Chen et al. (2017) is based on direct numerical simulation (DNS), i.e. solving the governing equations for the exact reconstructed geometry of the electrode. DNS by its very nature requires enormous computing resources and therefore both these studies only simulated a small fraction of an actual RFB electrode that while providing some useful insights, it is not necessarily representative. Furthermore, even for a relatively small domain like those in these works, optimization or parametric studies become extremely time-consuming as the governing equations need to be solved repeatedly. Finally,  only studied the effect of porosity on the RFB performance. There are other important structural features such as pore size distribution, fiber alignment, etc. that could potentially improve RFB performance that are yet to be discovered. In this study, we first introduce a new modeling framework to study the RFB at pore-scale with minimal computing resources. This framework is based on pore networks which has been previously demonstrated to reduce the simulation time by 4 to 5 orders of magnitudes compared to that of the LBM. 37 We then demonstrate the developed model in the context of a hydrogen-bromine RFB and explore how different microstructural features affect the overall RFB performance.

Model Development
Problem formulation.-The present study focuses on the cathode half-cell of a hydrogen bromine RFB with interdigitated flow configuration operating under steady-state conditions. One advantage of modeling this system is that the anode side could be neglected in the simulations as hydrogen oxidation is facile and therefore contributes little to the overall performance of the flow battery. 38,39 It should be emphasized that the insights obtained with the model can be translated to any RFB chemistry, of which many are currently under consideration. Figure 1a illustrates the schematic of a hydrogen-bromine flow battery. At the anode side, hydrogen gas H 2 is oxidized to form protons H + which then travel through the membrane to the cathode side where they react with bromine Br 2 to form hydrobromic acid HBr under the following electrochemical reaction. Figure 1b shows the cross-section view of the cathode compartment, which is comprised of several repeated units of flow channels for distributing bromine and hydrobromic acid, and ribs for collecting the generated current. Due to symmetry, the computational domain was chosen to be a half channel of bromine, a full rib, and a half channel of hydrobromic acid, as shown in Figure 1c. Here, the computational domain is only the porous electrode, and the flow channels, rib, and the membrane define the boundary conditions. Note that the flow channels are interdigitated, meaning that the bromine concentration remains constant along the channel, allowing constant Dirichlet boundary condition for concentration of bromine at inlet.
Pore network modeling.-Briefly, pore network modeling is mapping the porous domain onto a network of pores, which represent the void space, interconnected with arbitrary throats. 40 While this approach leads to loss of information at pore level, it enables the study of much larger porous domains with limited computing power that otherwise would be impossible if a full 3D finite element/volume method (FEM/FVM) had been employed to solve the governing equations. 41 This advantage of the pore network approach has recently made it an attractive alternative for studying how the microstructure affects the macroscopic properties of porous materials, which has been demonstrated in a few recent studies. [42][43][44] The basic assumption behind pore network modeling is that each pore is a well-mixed entity, i.e. within the pores, intensive properties such as pressure and concentration vary only slightly, if at all. This assumption is reasonably valid under a wide range of flow and transport regimes, including those in typical flow battery electrodes, and have been demonstrated in numerous studies. A recent study by Yang et al. is a good example in which the validity of this assumption is verified by comparison to experimental data and two independent direct numerical simulations using FEM and LBM. 45 Network generation.-A cubic pore network model was employed to represent the flow battery electrode. The network was assumed to be fully-connected, i.e. the coordination number of the internal pores is 6, and is equally spaced in x, y, and z, directions with a uniform centerto-center spacing of about 65 microns. A log-normal bimodal pore size distribution (PSD) based on the work of Zenyuk et al. (2016) was used, which reasonably represents the porous structure of SGL 25AA, a commercial gas diffusion layers that is widely used in flow battery testing. 46 Figure 2a shows a representative pore network model used in this study. Pores and throats are signified by spheres and cylinders, respectively, and the color legend shows the pore size distribution in the network. Figure 2b shows the histogram plots representing the pore/throat size distributions. Note that for the subsequent sensitivity analyses on structural features such as porosity, the structure of SGL 25AA was used as the base network and then modified accordingly to accommodate the requirements of that specific study.
Mathematical modeling.-The physics involved in the flow battery electrode at the cathode side are advection-diffusion of bromine, conduction of protons/electrons, and reaction of bromine and protons at the internal surface of the electrode. We also assume that the electrolyte is dilute enough and therefore, species transport does not affect the flow field. For simplicity, we also assume that the reaction is only taking place inside the pores and not throughout the connecting throats. This assumption allows for easier calculation of the flux between two adjacent pores. The introduced error due to underestimation of the active reaction area can be circumvented by adding half the surface area of the throat to each of the adjacent pores.
where ρ is the density of the electrolyte, N i is the number of neighbour pores to pore i, and finally u i j and S i j are the fluid velocity from pore i to j, and the cross-section area of the connecting throat, respectively. Assuming incompressible laminar flow, Hagen-Poiseuille equation can be used to express u i j based on pore pressures as defined by Where p i and p j are the respective pressures at pore i and j, and α i j = S i j /8πμ i j is the hydraulic conductance of the connecting throat of length i j . As for the boundary conditions, inlet flow rate, Q in , and discharge pressure p out = 0 were enforced at the inlet and outlet boundaries, respectively. No-flux condition was applied to the rest of the boundaries.
Bromine transport.-Writing the conservation of species for bromine around pore i gives where R i is the net reaction rate of bromine in pore i, κ i is the apparent rate constant, and m i j is the mass flux of bromine flowing from pore i to j which comprises both advective and diffusive terms. Normally in pore network models, the advective flux is discretized based on the first order upwind differencing scheme. Recently, we have demonstrated that this approach could be erroneous and developed a new pore-scale model based on the exact solution of 1D advection-diffusion in throats. 47 Based on their approach, the mass flux flowing from pore i to j is expressed by the following equation where c i and c j are concentrations of bromine at pore i and j, respectively, and Pe = u i j i j /D is the local Peclet number at the connecting throat. For boundary conditions, the concentration of bromine at the inlet was held constant at c in , while no-flux condition was applied to the rest of the boundaries. Note that the reaction rate of bromine is coupled with generated current within the electrode and therefore, the apparent rate constant is not truly a constant. Once proton/electron transport is explained, a relationship for calculating κ i based on relevant variables will be derived.
Proton/Electron transport.-Finally, writing a charge balance around pore i gives the governing equation for transport of protons as [6] where I i j is the charge flux flowing from pore i to j, and R p i is the net reaction rate of protons in pore i. The charge flux I i j is linearly proportional to the potential difference between the two pores, and therefore can be replaced with where φ i and φ j are the electric potential of the electrolyte at pore i and j, respectively, and β i j = σ / i j is the electrical conductance of the connecting throat where σ is the bulk electrical conductivity of the electrolyte, and is a weak function of concentration for strong acids such as HBr. The charge generation rate R p i was calculated using Butler-Volmer kinetics defined by where j 0 is the exchange current density, α c and α a are transfer coefficients for cathode and anode, respectively, where α a + α c = 1, A i is the internal surface area of pore i, c i is the concentration of bromine at pore i, c 0 is the reference concentration of bromine at which the open-circuit voltage V oc was measured, z is the number of electrons involved in the cathodic reaction, R is the universal gas constant, and T is the operating temperature. η c is the overpotential at the cathode side and is defined as the potential difference between the solid phase and the liquid phase: As explained previously, we assumed that the anode side contributes little activation polarization and IR loss. Accordingly, the voltage at the membrane should be roughly around φ = 0 V . However, since the electric conductivity of the membrane is on the order of O(10) S.m −1 , which is relatively low, the voltage loss across the membrane might not be negligible, specially when the electrode is operating far from equilibrium, i.e. V cell V oc . To resolve this issue, we considered an average voltage loss of φ m across the membrane and therefore, the boundary condition at the membrane interface becomes φ = φ m . For calculating φ m , one could write the Ohm's law across the membrane as defined by where R m is the ohmic resistance, and I m is the net current passing through the membrane. Since the electrode is operating at steady state Table I. Summary of the variables used in this study.

Variable Units Description
Current generation at pore i conditions, all protons passing through the membrane must get consumed somewhere within the electrode. Therefore, I m can be calculated by summing the proton consumption rate R p i over all internal pores: [11] Finally, the ohmic resistance of the membrane R m can be related to its conductance σ m through the Pouillet's law as defined by [12] where δ m and A m are the thickness and the cross-section area of the membrane, respectively. Note that the electrolyte voltage boundary condition at the membrane interface depends of the value of I m , and yet I m cannot be determined until after the simulation is finished. Hence, an iterative scheme is required, which is explained in section 0. For the rest of the boundaries, no-flux condition was applied. As for electron transport, since the electric conductance of the solid phase is on the order of O(10 3 ) S.m −1 , which is relatively high, we did not solve the electron conservation equation to lower the computation time even further. The electric potential of the solid phase was thus assumed to be constant throughout the electrode and equal to the applied voltage to the rib φ s = V cell . Finally, for calculating the apparent rate constant κ i , note that the proton consumption rate R p i is linearly proportional to the reaction rate of bromine R i through the Faraday constant F as described by the equation Combining Eqs. 4, 8, and 13, one could simply solve for the apparent rate constant κ i to obtain The variables and parameters used in this study and their respective units, values, and descriptions are summarized in Tables I and II, re-  Numerical solution.-The fluid flow in the electrode is governed by Eq. 2, which describes how pressure changes for each individual pore and its neighbors. Written for every pore in the network, the resulting system of linear equations can be readily solved to obtain the pressure field. Note that since the electrolyte was assumed to be dilute, fluid flow is not coupled with other transport mechanisms and therefore, can be solved separately. As for the bromine and proton transport, writing Eqs. 4 and 6 for every pore in the network results in two nonlinear systems of equations that need to be solved to obtain the concentration of bromine and the electric potential of the electrolyte, respectively. Note that since current generation is limited by reaction of bromine, the two systems of equations are coupled. An iterative scheme based on separately yet consecutively solving the two systems of equations was employed to handle the coupling between the two physics, in which the overpotential η c is initially guessed and then refined during each iteration. A detailed flowchart of this algorithm is illustrated in Figure 3. In the flowchart, the blue circles at the bottom right of each box refers to the corresponding equation for that step.
It is worth mentioning that the charge transport equation is highly nonlinear due to the exponential dependence of the current generation term on overpotential, and therefore the proposed explicit algorithm diverges quickly if the initial guess is not close to the correct solution. A relaxation factor ω was used to mitigate this issue such that φ k+1 where k is the iteration number. Furthermore, the coupling between the mass transport and charge transport equations is relatively weak when the cell is operating at close to the open-circuit voltage. Under such circumstances, the kinetic driving force η c is relatively small, leading to a low reaction rate and consequently an almost uniform concentration of bromine within the electrode. Hence, the two physics are weakly coupled and consequently, even with a relatively high relaxation factor, i.e. close to 1, convergence is guaranteed.
However, as the applied cell voltage deviates further from its opencircuit value, the coupling between the two physics becomes progressively stronger and a smaller relaxation factor is required to guarantee convergence. While the optimal value for the relaxation factor is a matter of trial and error and depends on several parameters such as the exchange current density, size of the electrode, etc., for the special

Applied voltage (V c )
Relaxation factor (ω) case that was studied in this paper, we recommend using the values reported in Table III.
Inferring fiber size in cubic networks.-Flow battery electrodes are typically made of fibrous materials. In this work, regular pore networks were used to model the porous domain, where pores are placed at the vertices of a regular cubic lattice. Therefore, since such networks are constructed using information about the void space, i.e. pores/throats, information related to the solid phase such as the average fiber diameter must be somehow inferred. Figure 4a shows a 2D cross-section of a sample cubic network, where pores, throats, and the solid phase, i.e. fibers, are shown in white, light blue, and black. Based on this figure, the cross-section of each individual fiber was defined to be the region colored in black enclosed by the pores on each of the square lattices. The total cross-sectional area of the fibers, A f , was then calculated using an image-based approach assuming spherical pores and cylindrical throats. Finally, assuming cylindrical fibers, their approximate average diameter was calculated using the following formula [15] where N f is the total number of fibers (as defined above) visible in the 2D cross-section of the network. Note that depending of which crosssection of the network is used, the calculated value for average fiber diameter slightly changes, but since the pore size was randomly distributed, the changes are insignificant across different cross-sections. When studying the effect of porosity, one needs to keep the fiber diameter relatively constant since the varying porosity is usually a result of physical compression, hence fiber diameter is not affected. Varying the porosity while maintaining a relatively constant fiber diameter is not trivial in cubic networks. Figure 4 is a demonstration of how this can be achieved, where a 2D cross-section of a cubic network is undergone a series of topological manipulation to achieve different microstructural properties. Here, d p is the average pore di-ameter, d f the average fiber diameter, δ the network spacing, and φ the porosity of the network. Note that for fair comparison, the macroscopic dimensions of all the five networks in this figure were kept constant at 400 μm by 400 μm. Based on Fig. 4b, increasing network spacing while keeping other parameters constant, leads to an increase in fiber diameter and a decrease in porosity. This is because when the spacing is increased, fewer pores can occupy the same area, which leads to a decrease in porosity and consequently an increase in the average fiber diameter. On the other hand, increasing the pore size while keeping other parameters constant, increases the porosity and therefore, decreases the average fiber diameter as shown in Figure 4c. Therefore, by carefully increasing both the spacing and the pore size, one could achieve different porosities while maintaining a fixed fiber diameter. This trial-and-error process is demonstrated in 4d and 4e, where two different porosities of φ = 59% and 49% were achieved at a fixed fiber diameter of d f = 36 μm.

Results and Discussion
Rather than performing a specific case study, the scope of this study was to introduce the use of pore network models for studying microstructural effects on the performance of RFBs. Nevertheless, we show that with no forced calibration, the results from our model agree very well with experimental data in the literature. As an example, we compared the results from our model with those extracted from an experimental study on a high-performance hydrogen bromine flow battery by Cho et al. (2012) under the same parameters, as shown in Figure 5a.
The only calibration applied was tuning the effective conductivity of the membrane as it was not reported by the reference article. Note that a volume-averaged continuum model "could" have achieved the same accuracy except in that case, one needs to either experimentally measured the effective properties such as tortuosity, permeability, conductivity, and active surface area or use them as fitting parameters to be able to match the results of an actual experiment. On the other hand, a pore network model does not rely on effective properties since it resolves the physics at pore scale. In fact, such properties can often be estimated from a pore network simulation, which can later be used as input parameters in a volume-averaged continuum model. Furthermore, the alternative method to capture the physics at pore scale is to run direct numerical simulations on 3D tomography images. If on average every pore in the tomography data is represented by 10 to 20 voxels in each direction, such simulations take roughly 3 to 4 orders of magnitude longer to finish. Besides, 3D tomography images are not always readily available. Conversely, pore network models are much quicker and only require an estimation of the pore size distribution of the material. We concede that these advantages come at the cost of accuracy. Nonetheless, it enables us to study much larger domains in a reasonable timeframe. Therefore, we believe that the advantages of pore network models over direct numerical simulations outweigh their relatively small loss in accuracy.
Finally, one might object that comparing macroscopic features such as the polarization curve is not sufficient for validating a pore-scale simulation, and that ideally, pore-scale concentration profiles should also be compared to fully validate the model. While this objection is valid, obtaining pore-scale concentration profiles from lab experiments is very difficult -if not impossible -due to the extremely small thickness of RFB electrodes, which are typically only a few hundred microns thick. In principle, it is possible to compare to direct numerical simulation (DNS) on tomography images, but the present model is based on a cubic geometry, so a direct comparison would not be helpful.
The results of this study are presented in the format of a series of sensitivity analyses on how the microstructural features of the flow battery electrode namely porosity (at constant fiber diameter) and fiber alignment affect its overall performance. Typical results obtained from the simulations are shown in Figure 5a and Figure 5b in form of polarization and power curves and plane-averaged heatmaps for concentration of bromine and electrolyte voltage. The concentration and voltage heatmaps shown in Figure 5b refer to an electrode with a porosity of ≈ 50%, plotted at four different operating voltages V c between 0.9 and 0 V with an exchange current density of j 0 = 0.5 A.m −2 . In the following subsections, the effect of porosity and fiber alignment on the overall performance of the electrode are presented and discussed in detail. Figure 6 shows how porosity affects electrode polarization and its power density. For a fair comparison, the electrode dimensions were held constant. Also note that varying porosity often is a result of different levels of compression and there- fore, the fiber diameter should also remain constant. To accommodate these requirements, different networks were constructed by the procedure explained in Sec. 0 by varying the network spacing, the number of pores in each direction, and finally the pore size distribution parameters, such that neither the electrode dimensions nor the fiber diameter change. Table IV provides a summary of the networks constructed at each porosity, where N i is the number of pores along the i-axis.

Effect of porosity.-
According to Figure 6a, as the porosity increases, more current is drawn from the electrode at a fixed voltage and therefore, the polarization curve shifts toward higher current densities. One could attribute this effect to higher reaction surface area of the electrode with higher porosity. However, this is not true since the networks were constructed such that the fiber diameter remains constant, which means varying porosity corresponds to different levels of compression of the same material. Therefore, since the electrode dimensions were also held constant, the electrode with lower porosity consists of more fibers and consequently has higher reaction surface area. Nevertheless, the electrode with lower porosity has a considerably lower permeability, which negatively affects the mass transfer, causing the polarization curve to rapidly plateau to a limiting current density. This effect can be observed in Figure 6a, where at relatively high porosities, the polarization curve does not plateau to a limiting current density, indicating a superior mass transfer compared to that at low porosities. The reason   that the negative effect of lower permeability at low porosity overwhelms the positive effect of higher reaction surface area is that flow rate is quadratically proportional to the pore diameter. Therefore, a 10% decrease in porosity, while increasing the reaction surface area by 10%, decreases the flow rate by 19%. Figure 6b shows a similar trend, where the electrode with higher porosity shows superior performance in terms of power density, compared to those with lower porosity.
A subtle, yet interesting feature of Figure 6b is that as porosity increases, the locus of the maximum power density tends to shift toward higher overpotentials. Figure 7 further illustrates this feature and shows that the maximum power density occurs at higher operating voltage for electrodes of higher porosity. According to this figure, as the porosity increases, the electrode must operate at higher overpotentials to reach the peak power density. Note that the kinetic and concentration driving forces balance out each other such that an increase in the former indirectly decreases the latter. For instance, as the kinetic driving force increases, the concentration driving force starts to decrease because of the increased reaction rate. However, at relatively high porosities, the improved mass transfer enables the electrode to operate at relatively higher overpotentials before the effects of mass transfer limitations show. Therefore, as the porosity increases, the locus of the maximum power density shifts toward higher overpotentials, i.e. cell voltages farther from V oc . Nevertheless, while the electrode with higher porosity provides higher power density, it must operate at a higher overpotential to reach its peak point, leading to a worse voltage efficiency. Hence, the higher power density of highly porous electrodes could be misleading. In fact, flow batteries are recommended to operate at only a fraction of their peak power density to minimize such losses. Nevertheless, at a constant flow rate, the electrode with higher porosity would reduce the pumping costs because of its higher permeability.
Note that the curve plateaus at relatively high porosities. This suggests that the overpotential where peak power occurs becomes independent of porosity. This is because at relatively high porosities, mass transfer is no longer the limiting step and therefore the distribution of  bromine is determined purely by kinetics rather than the hydrodynamic features of the electrode, i.e. porosity. Hence, further increasing the porosity does not significantly impact the locus of maximum power density. Finally, it should be noted that these results are based on the hydrogen bromine RFB system, which undergoes fast kinetics. For those with sluggish kinetics such as the all-vanadium RFB system, the reduction in active surface area per unit volume due to an increase in porosity and the subsequent performance loss could overshadow the gain due to the increased permeability. We hesitate to definitively make such conclusions, but still would like to emphasize that the trade-off between mass transfer and kinetics in such systems is more delicate and needs to be separately investigated.
Effect of fiber alignment.-Another microstructural feature of porous electrodes is directional anisotropy. For this study, we are interested specifically in capturing the effect of fiber alignment in different directions on the overall performance of the electrode. To incorporate the fiber alignment in the pore network model, a method developed by Gostick et al. (2007) was used in which pores with similar size are correlated and grouped along the desired direction. 48 Their method involves a tuning parameter β to control the degree of alignment and is defined as the number of adjacent pores affected by the correlation. A correlation distance of β = 0 indicates no alignment and therefore, represents a random fibrous electrode. We compared four different electrodes with different degrees of alignment as determined by the correlation distance parameter β. The perspective and top views of these electrodes are shown in Figure 8.
Since the problem in inherently symmetric with respect to the xaxis, we only studied the effect of alignment along the y-axis, i.e. perpendicular to the flow channels. We then evaluated their performance in terms of power density and plotted the results in Figure 9.
Based on this figure, fiber alignment has improved the maximum power density by a considerable amount of 33%. The increase in power density can be attributed to the fact that the aligned electrode has higher permeability along the alignment axis and therefore, has higher mass transfer coefficient due to the resulting higher flow rate under a fixed pressure drop. This effect has been previously verified experimentally by Kok et al. where they observed that electrospun electrodes produced with different degrees of anisotropy, by using higher collector drum rotation speeds, show increased permeability and diffusivity along the aligned axis. 49 Figure 10 further illustrates the effect of fiber alignment by comparing the concentration contours of bromine at different degrees of alignment. The red line in this figure indicates the longitudinal position at which the average concentration of bromine is 90% of its value at the inlet.
Based on this figure, fiber alignment causes the location of the red line to shift toward the outlet. Therefore, in the electrode with aligned fibers, more pores are operating at concentrations closer to that at the inlet, leading to a higher reaction rate and thus a higher power density. Although alignment appears to help the situation, it can be seen in Figure 9 that beyond a correlation distance β = 2 further aligning the electrode fibers (β > 2) does not further boost its performance. This is rather confusing since the extra fiber alignment should expectedly increase the permeability, improve the advective mass transport and consequently, the electrode should perform much better. To explain this observation, consider the permeability of the electrode at different degrees of fiber alignment shown in Figure 11. As can be seen, increasing the correlation distance beyond β = 2 does not significantly increase the permeability, and explaining the diminishing performance gains of the electrodes with a greater correlation distance.

Conclusions
In this work, a pore-scale model was developed based on pore networks, to study the multiphysics transport in redox flow batteries. To the best of our knowledge, this is the first time that pore network methodology is applied to study flow battery performance, though pore-network characterization of structural parameters has recently been reported. 30,31 The use of pore network approach enabled us to study a large domain spanning the entire thickness of the electrode, that is otherwise infeasible using direct numerical simulations such as  finite element or lattice Boltzmann methods. Therefore, the proposed framework is an ideal candidate to perform geometry optimizations at pore scale in a reasonable time. The model was demonstrated in the context of a hydrogen-bromine flow battery and was validated against the experimental data with the pore size distribution parameters as the only calibration tool but can be readily applied to any RFB chemistry. The effect of porosity at constant fiber diameter, and fiber alignment on electrode performance was studied. As for porosity, it was found that increasing the porosity at constant fiber diameter decreases the reaction surface area but increases permeability more significantly and therefore the overall performance monotonically increases. Having said that, it was also shown that as porosity increases, the loci of maximum power density is shifted toward higher overpotentials. Thus, the higher performance comes at the cost of lower voltage efficiency, which might not necessarily translate to a more economically feasible electrode. Finally, it was shown that electrodes with fibers aligned along the flow direction have better performance. However, the performance gain was shown to experience diminishing returns as the degree of alignment was increased. This was explained by the slower increase in permeability as a function of fiber alignment, which further emphasizes the importance of permeability among other key parameters. The effect of other microstructural features such as pore size distribution and inclusion of rigorous efficiency calculations including pumping costs is highly suggested to be studied as an extension to this work.