Theory of Multicomponent Phenomena in Cation-Exchange Membranes: Part III. Transport in Vanadium Redox-Flow-Battery Separators

In focusing our study on the transport driving forces in VRFB membranes (see Fig. 1 bottom), we consider the cell at quasi-steady-state in which the cell state of charge and the composition of the electrode electrolytes does not change significantly during a cycle. This work compliments previous VRFB modeling efforts that have focused on overall cell operation (see Fig. 1 top).^{10,12,13,25,26,34–40} This section first summarizes the VRFB system, including relevant electrochemical and chemical reactions. We then show how ohmic overpotentials and crossover of active species impact cell operation. The final sections formulate the transport and thermodynamic equations governing species flux and chemical potential drop across the membrane.

This study examines an all-vanadium redox-flow battery, pictured in Fig. 1. The system consists of electrolyte-filled porous electrodes that contain aqueous vanadium sulfate salts with sulfuric-acid supporting electrolyte (water, H₂O, protons, H⁺, and sulfate, SO₄²⁻, are abbreviated 0, H, and SO₄, respectively). The system is at a temperature of 295 K. The negative electrode, $\beta$, contains ${{\rm{V}}}^{2+}$ and ${{\rm{V}}}^{3+}$ (denoted V(II) and V(III), respectively, based on vanadium oxidation state) and the positive electrode, ${\epsilon },$, contains ${{\rm{VO}}}^{2+}$ and ${{\rm{VO}}}_{2}^{+}$ (denoted V(IV) and V(V) respectively). Vanadium in any oxidation state is denoted ${\rm{V}}({\rm{x}})$. As the cell discharges, the reactions are ${\rm{V}}\left({\rm{II}}\right)\to {\rm{V}}({\rm{III}})$ at the negative electrode and ${\rm{V}}\left({\rm{V}}\right)\to {\rm{V}}({\rm{IV}})$ at the positive electrode. Sulfate species associate with vanadium ions and with protons to form various ion-paired products.⁴¹ Table I details the electrochemical reactions and their half-cell potentials, the chemical reactions that occur when vanadium crosses over the membrane separating the electrodes (phase $\delta$), and the ion association equilibrium reactions.

Table I.  Electrochemical and chemical reaction and ion-pair equilibria in VRFBs.^22,41

Electrochemical Reactions
${{\rm{V}}}^{2+}\leftrightarrows {{\rm{V}}}^{3+}+{{\rm{e}}}^{-}$	${E}^{0}=-0.26V\,\mathrm{vs}\,\mathrm{SHE}$
${{\rm{VO}}}^{2+}+{{\rm{H}}}_{2}{\rm{O}}\leftrightarrows {{\rm{VO}}}_{2}^{+}+2{{\rm{H}}}^{+}+{{\rm{e}}}^{-}$	${E}^{0}=1.00V\,\mathrm{vs}\,\mathrm{SHE}$
Chemical Crossover Reactions
${{\rm{V}}}^{2+}+2{{\rm{VO}}}_{2}^{+}+2{{\rm{H}}}^{+}\to 3{{\rm{VO}}}^{2+}+{{\rm{H}}}_{2}{\rm{O}}$	At Positive Electrode
${{\rm{V}}}^{3+}+{{\rm{VO}}}_{2}^{+}\to 2{{\rm{VO}}}^{2+}$	At Positive Electrode
${{\rm{VO}}}^{2+}+{{\rm{V}}}^{2+}+2{{\rm{H}}}^{+}\to 2{{\rm{V}}}^{3+}+{{\rm{H}}}_{2}{\rm{O}}$	At Negative Electrode
${{\rm{VO}}}_{2}^{+}+2{{\rm{V}}}^{2+}+4{{\rm{H}}}^{+}\to 3{{\rm{V}}}^{3+}+2{{\rm{H}}}_{2}{\rm{O}}$	At Negative Electrode
Ion Association Equilibria
${{\rm{H}}}^{+}+{{\rm{SO}}}_{4}^{2-}\rightleftarrows {{\rm{HSO}}}_{4}^{-}$
${{\rm{HSO}}}_{4}^{-}+{{\rm{V}}}^{3+}\rightleftarrows {{\rm{HSO}}}_{4}{{\rm{V}}}^{2+}$
${{\rm{HSO}}}_{4}^{-}+{{\rm{VO}}}^{2+}\rightleftarrows {{\rm{HSO}}}_{4}{{\rm{VO}}}^{+}$

Cell performance is typically characterized by the power density on discharge, ${{\rm{\Psi }}}_{{\rm{d}}}\,$, and the round-trip energy efficiency, ${\varepsilon }_{{\rm{e}}}$.⁴² ${{\rm{\Psi }}}_{{\rm{d}}}$ is the product of the discharge current density ${i}_{{\rm{d}}}$ and the cell potential on discharge ${V}_{{\rm{d}}}$⁴³

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{d}}}={i}_{{\rm{d}}}{V}_{{\rm{d}}}\end{eqnarray} \tag{ 1 }$$

The energy efficiency is the ratio of the integrated cell power density on discharge to charge and is typically treated as the product of the voltage efficiency, ${\varepsilon }_{{\rm{v}}}$, and the Coulombic efficiency, ${\varepsilon }_{{\rm{q}}}$,⁶

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{e}}}=\left(\displaystyle \frac{{V}_{{\rm{d}}}}{{V}_{{\rm{c}}}}\right)\,\left(\displaystyle \frac{{\int }_{{\rm{d}}}{i}_{{\rm{d}}}dt}{{\int }_{{\rm{c}}}{i}_{{\rm{c}}}dt}\right)={\varepsilon }_{{\rm{v}}}{\varepsilon }_{{\rm{q}}}\end{eqnarray} \tag{ 2 }$$

where $t$ is time, the subscript ${\rm{c}}$ or ${\rm{d}}$ denotes charge or discharge, and the cell voltage is approximated as constant during charge or discharge. ${\varepsilon }_{{\rm{v}}}$ characterizes efficiency losses due to cell overpotentials (e.g. kinetic and ohmic) and decreases with increasing current density.³ ${\varepsilon }_{{\rm{q}}}$ quantifies efficiency losses due to vanadium passing through the membrane between electrode solutions, mixing, and reacting chemically rather than electrochemically. ${\varepsilon }_{{\rm{q}}}$ tends to increase with rising current density because the total charge extracted from the cell generally grows compared to the amount of charge lost to chemical reactions.³ A typical operating goal is maximizing power density at a target energy efficiency.⁴² The next two sections present the governing equations for cell potential that are required to calculate power density and voltage efficiency, and the state-of-charge loss during cycling that dictates Coulombic efficiency.

Cell Potential

The cell potential is the difference in the electrochemical potential of the electrons ${\mu }_{{{\rm{e}}}^{-}}$ in the metal at the negative electrode $\alpha ^{\prime}$ and the positive electrode $\alpha ^{\prime\prime}$; these are related to the electrochemical potential of species $i$, ${\mu }_{i}^{\phi }$, participating in the oxidation/reduction reactions in electrode $\phi$ as outlined in Table I⁴³

$$\begin{eqnarray}&&\begin{array}{l}FV={\mu }_{{{\rm{e}}}^{-}}^{\alpha ^{\prime} }-{\mu }_{{{\rm{e}}}^{-}}^{\alpha ^{\prime\prime} }=\left({\mu }_{{\rm{V}}\left(\mathrm{II}\right)}^{\beta }-{\mu }_{{\rm{V}}\left(\mathrm{III}\right)}^{\beta }-F{\eta }^{\beta }\right)\\ \,\,\,-\,\left({\mu }_{{\rm{V}}\left(\mathrm{IV}\right)}^{\epsilon }+{\mu }_{0}^{\epsilon }-{\mu }_{{\rm{V}}\left({\rm{V}}\right)}^{\epsilon }-2{\mu }_{{\rm{H}}}^{\epsilon }-F{\eta }^{\epsilon }\right)\end{array}\end{eqnarray} \tag{ 3 }$$

where $F$ is Faraday's constant and ${\eta }^{\phi }$ is ionic potential drop between the metal and the electrolyte in electrode $\phi$ due to mass-transport resistances, surface overpotentials, or other non-membrane cell resistances.

To put Eq. 1 in terms of the electrochemical potentials of species in the membrane, chemical equilibrium requires that the electrochemical potential of charged species (and the chemical potential of neutral species) must be equal at the interface of electrode phase $\phi$ and membrane phase $\delta$⁴⁴

$$\begin{eqnarray}&&{\mu }_{i}^{\delta }={\mu }_{i}^{\phi }\end{eqnarray} \tag{ 4 }$$

where ${\mu }_{i}$ is a function of chemical variables—temperature, pressure, composition—and, for charged species, the ionic potential.

It is convenient to group electrostatic dependences into a single variable with the remaining variables depending only on composition, temperature, and pressure. We reference the electrochemical potential of species $i$ to the electrochemical potential of the membrane ${\rm{M}}$ to define

$$\begin{eqnarray}&&{\mu }_{i,{\rm{M}}}^{\delta }={\mu }_{i}^{\delta }-\frac{{z}_{i}}{{z}_{{\rm{M}}}}{\mu }_{{\rm{M}}}^{\delta }\end{eqnarray} \tag{ 5 }$$

which is an electroneutral pairing of species and is independent of electrostatic potential.^43,45 To track the electrical state through the membrane, we set the ionic potential ${\rm{\Phi }}$ equal to the electrochemical potential of the membrane, ${z}_{{\rm{M}}}F{\rm{\Phi }}={\mu }_{{\rm{M}}}^{\delta }$. This assignment is equivalent to measuring the ionic potential using a hypothetical reference electrode selective to the ionic groups on the membrane at the same conditions.⁴³ Using Eq. 4 and ${\mu }_{i,{\rm{M}}}^{\delta }$, the cell potential is then

$$\begin{eqnarray}&&V=U+\eta +{\rm{\Delta }}{\rm{\Phi }}\end{eqnarray} \tag{ 6 }$$

where $\eta ={\eta }^{{\epsilon }}-{\eta }^{\beta }$ groups non-membrane overpotentials in the cell. $U$ is the equilibrium cell potential

$$\begin{eqnarray}&&\begin{array}{l}FU=\left[\left({\mu }_{{\rm{V}}\left(\mathrm{II}\right),{\rm{M}}}^{\delta ^{\prime} }-{\mu }_{{\rm{V}}\left(\mathrm{III}\right),{\rm{M}}}^{\delta ^{\prime} }\right)-\left({\mu }_{{\rm{V}}\left(\mathrm{IV}\right),{\rm{M}}}^{\delta ^{\prime\prime} }+{\mu }_{0}^{\delta ^{\prime\prime} }\right.\right.\\ \,-\,\left.\left.{\mu }_{{\rm{V}}\left({\rm{V}}\right),{\rm{M}}}^{\delta ^{\prime\prime} }-2{\mu }_{H,M}^{\delta ^{\prime\prime} }\right)\right]\end{array}\end{eqnarray} \tag{ 7 }$$

where $\delta ^{\prime}$ and $\delta ^{\prime\prime}$ are the membrane phase at the negative and positive electrode interfaces, respectively. ${\rm{\Delta }}{\rm{\Phi }}={{\rm{\Phi }}}^{\delta ^{\prime\prime} }-{{\rm{\Phi }}}^{\delta ^{\prime} }$ is the potential drop across the membrane due to ohmic and diffusion overpotentials. During discharge (charge), $\eta$ and ${\rm{\Delta }}{\rm{\Phi }}$ become more negative (positive) with increasing current density, causing the cell voltage efficiency to decrease.³ We set $U$ to the measured open-circuit potential (=1.42 V)³ by neglecting the diffusion potential, which tends to be small.⁴³ η is modeled as an ohmic potential drop with resistance ${R}_{{\rm{\Omega }}}$

$$\begin{eqnarray}&&\eta =\pm | i| {R}_{{\rm{\Omega }}}\end{eqnarray} \tag{ 8 }$$

where $\eta$ is negative for discharge and positive for charge.

State-of-charge loss

Vanadium crossover causes the cell to self-discharge. The cell state of charge, $SOC$, characterizes the extent the cell is in the charged state; it is the ratio of the sum of the moles of vanadium in the charged state (namely, ${\rm{V}}({\rm{II}})$ in the negative electrode and ${\rm{V}}({\rm{V}})$ in the positive electrode) to the total moles of vanadium in the system

$$\begin{eqnarray}&&SOC=\displaystyle \frac{{n}_{{\rm{V}}\left({\rm{II}}\right)}^{\beta }+{n}_{{\rm{V}}\left({\rm{V}}\right)}^{{\epsilon }}}{{n}_{{\rm{V}}\left({\rm{II}}\right)}^{\beta }+{n}_{{\rm{V}}\left({\rm{III}}\right)}^{\beta }+{n}_{{\rm{V}}\left({\rm{V}}\right)}^{{\epsilon }}+{n}_{{\rm{V}}\left({\rm{IV}}\right)}^{{\epsilon }}}\end{eqnarray} \tag{ 9 }$$

where ${n}_{i}^{\phi }$ is the moles of $i$ in electrode $\phi .$

The rate at which the state of charge changes during discharging or charging is a function of the cell current density $i$ and the crossover current density ${i}_{{\rm{x}}}$^1–3

$$\begin{eqnarray}&&\displaystyle \frac{F{n}_{{\rm{t}}}}{2A}\displaystyle \frac{dSOC}{dt}=-i-{i}_{{\rm{x}}}\end{eqnarray} \tag{ 10 }$$

where $A$ is the active area of the membrane and ${n}_{{\rm{t}}}$ is the total moles of vanadium in the system (i.e. ${n}_{{\rm{V}}\left({\rm{II}}\right)}^{\beta }+{n}_{{\rm{V}}\left({\rm{III}}\right)}^{\beta }+{n}_{{\rm{V}}\left({\rm{V}}\right)}^{{\epsilon }}+{n}_{{\rm{V}}\left({\rm{IV}}\right)}^{{\epsilon }}$). ${i}_{{\rm{x}}}$ characterizes how vanadium crossover decreases the state of charge of the cell and is a function of the flux of vanadium species $i$ through the membrane, ${N}_{i}$, according to^1–3

$$\begin{eqnarray}&&{i}_{{\rm{x}}}=F\displaystyle \frac{3{N}_{{\rm{V}}\left({\rm{II}}\right)}-3{N}_{{\rm{V}}\left({\rm{V}}\right)}+{N}_{{\rm{V}}\left({\rm{III}}\right)}-{N}_{{\rm{V}}\left({\rm{IV}}\right)}}{2}.\end{eqnarray} \tag{ 11 }$$

The membrane/negative electrode interface is set at $z=0$ and the membrane/positive electrode interface is at $z={l}_{{\rm{M}}}$ so that ${N}_{i}$ is positive for species that move towards the positive electrode (e.g. ${\rm{V}}\left({\rm{II}}\right)$ and ${\rm{V}}({\rm{III}})$) and negative for species that move toward the negative electrode (e.g. ${\rm{V}}({\rm{IV}})$ and ${\rm{V}}({\rm{V}})$). With this convention for the direction of flux, the ionic current density $i$ through the membrane is positive on discharge and is also the electronic current leaving the positive electrode normalized by $A$ (i.e. positive on discharge).¹ A subsequent section shows that ${N}_{i}$ is constant across the membrane.

The first term on the right side of Eq. 10 is the consumption/production of species due to electrochemical charge/discharge, the second term is due to the loss of ${\rm{V}}({\rm{II}})$ and ${\rm{V}}({\rm{V}})$ from crossover to the opposing electrode and from reactions with vanadium species that transport from the opposing electrode (Table I specifies the stoichiometry incorporated in Eq. 11). Based on the low concentration of absorbed vanadium in the membrane, the chemical reactions in Table I occur in the electrode solutions and not in the membrane. The supporting electrolyte provides an excess of protons, and the high concentration of vanadium in the electrode solution rapidly reacts with any vanadium that crosses the membrane into the electrode.²² The factor of 1/2 on the left side of Eq. 10 results from considering the combined state of charge of the positive and negative electrodes. ${i}_{{\rm{x}}}$ is a function of current density and cell state of charge (see Eq. 11).^13,36

As a result of vanadium crossover, the amount of charge extracted from the cell during discharge, ${\int }_{{\rm{d}}}{i}_{{\rm{d}}}Adt$, is less than the amount put into the cell to recharge it to its original state of charge during charging, $-{\int }_{{\rm{c}}}{i}_{{\rm{c}}}Adt$.¹ Integrating Eq. 10 over a discharge/charge cycle with no net change in $SOC$ for a constant discharge and constant charge current density and crossover current, ${i}_{{\rm{x}}}$, and solving for ${\varepsilon }_{{\rm{q}}}$ gives¹

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{q}}}=\displaystyle \frac{1+\displaystyle \frac{{i}_{{\rm{x}},{\rm{c}}}}{{i}_{{\rm{c}}}}}{1+\displaystyle \frac{{i}_{{\rm{x}},{\rm{d}}}}{{i}_{{\rm{d}}}}}\end{eqnarray} \tag{ 12 }$$

where ${i}_{{\rm{c}}}$ is negative and ${i}_{{\rm{d}}}$, ${i}_{{\rm{x}},{\rm{c}}}$, ${i}_{{\rm{x}},{\rm{d}}}$ are positive. Because ${i}_{{\rm{x}}}$ is a function cell state of charge, the condition for Eq. 12 of constant ${i}_{{\rm{x}}}$ is achieved by a small change to $SOC$ over charge/discharge.

The following two sections discuss how the electrolyte partitions into and how it transports across the membrane. These expressions provide for calculation of ${\varepsilon }_{{\rm{q}}}$, ${\varepsilon }_{{\rm{v}}}$, and ${\rm{\Delta }}{\rm{\Phi }}.$ Although the following expressions are general for any flow-battery membrane, the specific microscale theories we use to calculate thermodynamic and transport coefficient are specific to phase-separated polymer cation-exchange membrane that Parts I and II present.^32,33 The most common class of these materials are perfluorinated sulfonic-acid (PFSA) ionomers, such as Nafion.⁴⁶ These materials consist of nanoscale, interconnected hydrophilic domains filled with aqueous electrolyte solution and side-chains terminated with anionic moieties that are tethered to the PFSA polymer.⁴⁶ Surrounding these domains are hydrophobic polytetrafluoroethylene backbone that provides structural support.⁴⁶

Membrane thermodynamics

Requisite for species to transport across the membrane, they must sorb from the electrode electrolyte into the membrane as Eq. 4 specifies. Chemical thermodynamics quantifies how these potentials are related to measurable quantities. The chemical potential of a neutral species is⁴⁴

$$\begin{eqnarray}&&{\mu }_{i}={\mu }_{i}^{\theta }+RT\,\mathrm{ln}\,{x}_{i}+{\mu }_{i}^{{\rm{ex}}}\end{eqnarray} \tag{ 13 }$$

where $R$ is the universal gas constant, $T$ is temperature, and ${x}_{i}$ is the mole fraction of $i$ ($=\sum _{j}{n}_{i}/{n}_{j}$ where ${n}_{j}$ is the moles of species $j$). The mole fraction of charged polymer groups in the membrane is ${x}_{{\rm{M}}}$ and the charge for most cation-exchange membranes is ${z}_{{\rm{M}}}=-1$. The reference chemical potential of $i$, ${\mu }_{i}^{{\rm{\theta }}}$, is for a hypothetical ideal solution at unit mole fraction ${x}_{i}$ and is a function of temperature and pressure.⁴⁴ The second term on the right side of Eq. 13 is the ideal solution contribution.⁴⁴ In this paper, the ideal-solution contribution is specified a for fully-dissociated electrolyte (i.e. the Experimental Construct discussed in Part I).³² The excess chemical potential accounts for ion/ion, ion/membrane, solvent/membrane, and solvent/ion interactions including ion-pair formation (see Table I), electrostatic interactions, ion solvation, steric confinement, and membrane swelling. Part I details calculation of the excess chemical potential ${\mu }_{i}^{{\rm{ex}}}.$³²

For a charged species, the electrochemical potential has the same form as Eq. 13 but further depends on the electrical state of the phase and charge number ${z}_{i}.$⁴³ However, expressing the electrochemical chemical potential using only neutral pairings of two charged of species (i.e. ${\mu }_{i}-\displaystyle \frac{{z}_{i}}{{z}_{j}}{\mu }_{j}$) eliminates the dependence on electrical potential.⁴³ Substitution of these neutral pairings of Eq. 13 into Eq. 4 relates species mole fraction in the membrane and electrode electrolyte⁴⁷

$$\begin{eqnarray}&&\displaystyle \frac{{x}_{i}^{\delta }{\left({x}_{j}^{\delta }\right)}^{-\displaystyle \frac{{z}_{i}}{{z}_{j}}}}{{x}_{i}^{\phi }{\left({x}_{j}^{\phi }\right)}^{-\displaystyle \frac{{z}_{i}}{{z}_{j}}}}=\exp \left[\displaystyle \frac{\left({\mu }_{i}^{{\rm{ex}},\phi }-{\mu }_{i}^{{\rm{ex}},\delta }\right)-\displaystyle \frac{{z}_{i}}{{z}_{j}}\left({\mu }_{j}^{{\rm{ex}},\phi }-{\mu }_{j}^{{\rm{ex}},\delta }\right)}{RT}\right]={{\rm{\Gamma }}}_{ij}\,\end{eqnarray} \tag{ 14 }$$

For neutral species (e.g. water) ${z}_{i}=0$ and, therefore, the $j$th species does not need to be specified in Eq. 14. Electroneutrality inside the membrane adds a constraint that fully specifies ion partitioning,⁴⁴

$$\begin{eqnarray}&&F\displaystyle \displaystyle \sum _{i}{z}_{i}{x}_{i}^{\delta }=0\end{eqnarray} \tag{ 15 }$$

If oppositely charged ions are used for $i$ and $j$ (such as the anion and cation of a salt) and ${{\rm{\Gamma }}}_{ij}$ is unity, the oppositely charged pairs $i$ and $j$ partition into the membrane according to ideal Donnan equilibrium.⁴⁷ If ${{\rm{\Gamma }}}_{ij}\gt 1$, $i$ and $j$ favorably partition into the membrane whereas if partitioning is unfavorable, ${{\rm{\Gamma }}}_{ij}\lt 1.$⁴⁷ For context, Part I shows that for a Nafion cation-exchange membrane in an aqueous hydrobromic-acid solution at concentration of 5 mol kg⁻¹-solvent, the term on the right side of Eq. 14 is ∼0.75 for HBr and 0.95 for water.³²

Upon sorption into the membrane, species transport according to their molar flux ${N}_{i}$ across the membrane ($z$-direction in the 1-D model). The molar flux obeys species conservation at steady-state⁴⁸

$$\begin{eqnarray}&&-\displaystyle \frac{\partial {N}_{i}}{\partial z}=0\end{eqnarray} \tag{ 16 }$$

Although the concentration of species in the electrodes changes during cell operation, transport in the membrane is assumed to be at steady-state because the concentration and pressure gradients in the membrane develop much faster (i.e. pseudo steady-state approximation).⁴⁸Absent temperature and hydraulic-pressure gradients, which are typically small, transport is driven by electrochemical potential gradients. Friction between species introduces drag that reduces total flux. Stefan-Maxwell-Onsager theory describes the balance between driving force and drag on $i$ not equal to ${\rm{M}}$ as^43,49,50

$$\begin{eqnarray}&&{c}_{i}\displaystyle \frac{\partial {\mu }_{i}}{\partial z}=\displaystyle \displaystyle \sum _{j\ne i,{\rm{M}}}\displaystyle \frac{RT}{{{\mathfrak{D}}}_{ij}{c}_{{\rm{T}}}}({c}_{i}{N}_{j}-{c}_{j}{N}_{i})-\displaystyle \frac{RT{c}_{{\rm{M}}}}{{{\mathfrak{D}}}_{i{\rm{M}}}{c}_{{\rm{T}}}}{N}_{i},\end{eqnarray} \tag{ 17 }$$

and for the membrane as

$$\begin{eqnarray}&&{c}_{{\rm{M}}}\displaystyle \frac{\partial {\mu }_{{\rm{M}}}}{\partial z}=\displaystyle \displaystyle \sum _{j\ne {\rm{M}}}\displaystyle \frac{RT{c}_{{\rm{M}}}}{{{\mathfrak{D}}}_{j{\rm{M}}}{c}_{{\rm{T}}}}{N}_{j}\end{eqnarray} \tag{ 18 }$$

where ${c}_{i}$ is the molar concentration of species $i$ and ${c}_{{\rm{T}}}$ is the total molar concentration in the membrane. (For convenience in this section, we do not superscript variables for quantities in the membrane phase $\delta$). As discussed in the preceding section, ${\mu }_{i}$ is a function of temperature, pressure, composition, and ionic potential ${\rm{\Phi }}$. ${{\mathfrak{D}}}_{ij}$ is the binary diffusion coefficient between species $i$ and $j$. ${{\mathfrak{D}}}_{j{\rm{M}}}$ is related to the friction coefficient between species $i$ and the membrane, ${K}_{i{\rm{M}}}$, that is described in Part II, according to ${{\mathfrak{D}}}_{i{\rm{M}}}={c}_{i}{c}_{{\rm{M}}}RT/{c}_{{\rm{T}}}{K}_{i{\rm{M}}}.$³³ The membrane is stationary (i.e. is the reference velocity) and the membrane concentration ${c}_{{\rm{M}}}$ is set as the molar concentration of charged groups on the polymer. Onsager reciprocal relations dictate that the diffusion coefficients are symmetric so that ${{\mathfrak{D}}}_{ij}={{\mathfrak{D}}}_{ji}$. Consequently, for a VRFB with eight species (water, sulfate ions, protons, the membrane, and four vanadium species), there are 28 transport coefficients. In this paper, Eqs. 17 and 18 consider ionic species as fully dissociated (i.e. uses the Experimental Construct discussed in Part II).³³ Transport coefficients ${{\mathfrak{D}}}_{ij}$ include the effect of ion-pair formation to ensure that this model is consistent with the various ion-paired species that exist. Part II outlines calculation of ${{\mathfrak{D}}}_{ij}$ and ${K}_{i{\rm{M}}}$ as a function of membrane water content and ion concentration.³³

Equations 17 and 18 are rigorous but inconvenient because they involve gradients in electrochemical potential that are not readily characterized and frame the driving force in terms of species fluxes, whereas experiments measure fluxes under applied forces. Appendix A shows that for constant total molar concentration, ${c}_{{\rm{T}}},$ eliminating the electrostatic dependence of the driving forces (i.e. making the substituting ${\mu }_{i}={\mu }_{i,{\rm{M}}}+\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}}{\mu }_{{\rm{M}}}$ as shown in Eq. 5),⁴⁵ and expanding chemical potential gradients in terms of composition variables gives the flux of species $i$ as the sum of migration and diffusion terms:

$$\begin{eqnarray}&&{N}_{i}=\displaystyle \frac{{t}_{i}i}{{z}_{i}F}-\displaystyle \displaystyle \sum _{j\ne {\rm{M}},n}{D}_{ij}\displaystyle \frac{\partial {c}_{j,{\rm{M}}}}{\partial z}\end{eqnarray} \tag{ 19 }$$

for $i\ne {\rm{M}},n,$ where $n$ is a reference species with charge number different from ${\rm{M}},$ ${t}_{i}$ is the transference number of species $i,$ and ${D}_{ij}$ are multicomponent diffusion coefficient between species $i$ and $j$ based on concentration driving forces. We assign protons as reference species $n.$ Although ${t}_{0}$ and ${z}_{0}$ are zero for water, the ratio ${t}_{0}/{z}_{0}$ is definite and equal to the electroosmotic coefficient $\xi .$ $\partial {c}_{j,{\rm{M}}}/\partial z$ is the concentration gradient of a (potentially hypothetical) salt consisting of $| {z}_{{\rm{M}}}|$ ions $j$ and $| {z}_{j}|$ charged polymer groups ${\rm{M}}$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{j,{\rm{M}}}}{\partial z}=\displaystyle \frac{1}{{s}_{j,{\rm{M}}}}\left(\displaystyle \frac{\partial {c}_{j}}{\partial z}-\displaystyle \frac{{z}_{j}{c}_{j}}{{z}_{{\rm{M}}}{c}_{{\rm{M}}}}\displaystyle \frac{\partial {c}_{{\rm{M}}}}{\partial z}\right)\end{eqnarray} \tag{ 20 }$$

where ${s}_{j,{\rm{M}}}$ is the number of particles constituting $j,$ ${\rm{M}}$ (${s}_{j,{\rm{M}}}=\left|{z}_{j}\right|+| {z}_{{\rm{M}}}|$ for ${z}_{j}\ne 0$). For neutral species (e.g. water), ${s}_{j,{\rm{M}}}=1$ so that the concentration gradient of $j,$ ${\rm{M}}$ is simply the concentration gradient of $j$ (i.e. $\partial {c}_{j,{\rm{M}}}/\partial z=\partial {c}_{j}/\partial z$). This formalism generalizes the treatment of concentration gradients in binary electrolytes⁴³ to an arbitrary number of species; the concentration gradient of $j,$ ${\rm{M}}$ reduces to the concentration gradient of the salt for a binary electrolyte with $\left|{z}_{-}\right|$ cations and $| {z}_{+}|$ anions.

Appendix A shows how ${D}_{ij}$ is related to the binary diffusion coefficients ${{\mathfrak{D}}}_{ij}$'s and chemical potential gradients. Specifically, ${D}_{ij}$'s are components of the $N-2$ by $N-2\,$diffusion matrix ${\boldsymbol{D}}$ that is the product of the matrix inverse of ${\boldsymbol{B}},$ an $N-2$ by $N-2$ matrix containing transport coefficients, and ${\boldsymbol{\chi }},$ an $N-2$ by $N-2$ matrix containing multicomponent thermodynamic factors

$$\begin{eqnarray}&&{\boldsymbol{D}}={{\boldsymbol{B}}}^{-1}{\boldsymbol{\chi }}{\boldsymbol{S}}\end{eqnarray} \tag{ 21 }$$

where ${\boldsymbol{S}}$ is a diagonal matrix with entries ${s}_{j,{\rm{M}}},$ and the entries of ${\boldsymbol{B}}$ for $i\ne j$ are

$$\begin{eqnarray}&&{B}_{ij}={x}_{i}\left[\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}{{\mathfrak{D}}}_{j{\rm{M}}}}+\displaystyle \frac{{z}_{j}}{{z}_{n}}\left(\displaystyle \frac{1}{{{\mathfrak{D}}}_{in}}-\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}{{\mathfrak{D}}}_{n{\rm{M}}}}\right)-\displaystyle \frac{1}{{{\mathfrak{D}}}_{ij}}\right]\end{eqnarray} \tag{ 22 }$$

and for $i=j,$

$$\begin{eqnarray}&&{B}_{ii}=\displaystyle \frac{{z}_{i}{x}_{i}}{{z}_{n}}\left(\displaystyle \frac{{z}_{n}}{{z}_{{\rm{M}}}{{\mathfrak{D}}}_{i{\rm{M}}}}+\displaystyle \frac{1}{{{\mathfrak{D}}}_{in}}-\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}{{\mathfrak{D}}}_{n{\rm{M}}}}\right)+\displaystyle \displaystyle \sum _{k\ne i}\displaystyle \frac{{x}_{k}}{{{\mathfrak{D}}}_{ik}}.\end{eqnarray} \tag{ 23 }$$

${{\boldsymbol{B}}}^{-1}$ is a matrix of thermodynamic diffusion coefficients. The entries of ${\boldsymbol{\chi }}$ are

$$\begin{eqnarray}&&{{\rm{\chi }}}_{ij}={\delta }_{ij}+\displaystyle \frac{{x}_{i}}{{x}_{j}}\displaystyle \frac{\partial {\mu }_{i,{\rm{M}}}^{{\rm{ex}}}}{\partial RT\,\mathrm{ln}({x}_{j}/{x}_{{\rm{M}}}^{{z}_{j}/{z}_{{\rm{M}}}})}\end{eqnarray} \tag{ 24 }$$

where ${\delta }_{ij}$ is the Kronecker delta function (${\delta }_{ii}=1$ and ${\delta }_{ij}=0$ for $i\ne j$). Equation 21 reduces to a well-established form for the salt diffusion coefficients in the case of a binary electrolyte.⁴³

The $N-2$ vector ${\boldsymbol{t}}$ of charge number-normalized transference numbers ${t}_{i}/{z}_{i}$ is

$$\begin{eqnarray}&&{\boldsymbol{t}}={{\boldsymbol{B}}}^{-1}{\boldsymbol{\varsigma }}\end{eqnarray} \tag{ 25 }$$

where ${\boldsymbol{\varsigma }}$ is an $N-2$ vector with entries

$$\begin{eqnarray}&&{\varsigma }_{i}=\displaystyle \frac{{x}_{i}}{{z}_{n}}\left(\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}{{\mathfrak{D}}}_{n{\rm{M}}}}-\displaystyle \frac{1}{{{\mathfrak{D}}}_{in}}\right).\end{eqnarray} \tag{ 26 }$$

In this approach, ${t}_{i}$ and ${D}_{ij}$ use a reference velocity of the membrane (i.e. ${N}_{{\rm{M}}}=0$), unlike other research that uses water velocity as the frame of reference, see Ref. 1. Due to this choice of the laboratory frame of reference, Eq. 19 should not include additional terms for solvent convection. The choice for species $n$ does not change the flux ${N}_{i}$ but does change the value of ${D}_{ij}.$ Since Eq. 19 is applicable to $i\ne n,$ we specify the flux of $n$ using the current density and fluxes of species $i\ne n$ (i.e. ${N}_{n}=(i/F-\displaystyle \sum _{i\ne {\rm{M}},n}{N}_{i}{z}_{i})/{z}_{n}$).

For negligible mass-transport resistance between the electrode electrolyte and membrane, chemical equilibria (i.e. Eqs. 14 and 15) specifies ${c}_{i,{\rm{M}}}$ at the interface of the membrane with the two electrodes. Solving Eqs. 16 and 19 using these concentrations as boundary conditions gives the concentration profiles and fluxes of species across the membrane for a set current density.

With specified fluxes, Eq. 18 calculates the ionic potential gradient because it is proportional to the electrochemical potential of the membrane,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}=\displaystyle \frac{RT}{{z}_{{\rm{M}}}F{c}_{{\rm{T}}}}\displaystyle \sum _{j\ne {\rm{M}}}\displaystyle \frac{{N}_{j}}{{{\mathfrak{D}}}_{j{\rm{M}}}}.\end{eqnarray} \tag{ 27 }$$

Integration of Eq. 27 across the membrane shows that ${\rm{\Delta }}{\rm{\Phi }}$ is a function of applied current density (ohmic overpotential) and flux of species at the open-circuit cell potential (diffusion overpotential). In the absence of concentration gradients, substitution of Eq. 19 into Eq. 27 identifies the membrane conductivity, $\kappa ,$ as

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{\partial {\rm{\Phi }}}{\partial z}=-i{\kappa }^{-1}=-i\left\{-\displaystyle \frac{RT}{{c}_{{\rm{T}}}{F}^{2}{z}_{{\rm{M}}}}\left[\displaystyle \frac{1}{{z}_{n}{{\mathfrak{D}}}_{n{\rm{M}}}}\right.\right.\\ \,\,\,\,+\,\left.\left.\displaystyle \displaystyle \sum _{j\ne {\rm{M}},n}\displaystyle \frac{{t}_{j}}{{z}_{j}}\left(\displaystyle \frac{1}{{{\mathfrak{D}}}_{j{\rm{M}}}}-\displaystyle \frac{{z}_{j}}{{z}_{n}{{\mathfrak{D}}}_{n{\rm{M}}}}\right)\right\}\right].\end{array}\end{eqnarray} \tag{ 28 }$$

As $\kappa$ increases for a fixed current density and concentration gradient, the potential drop across the membrane decreases.

Parts I and II discuss calculation of ${\mu }_{i}^{{\rm{ex}}}$ and ${{\mathfrak{D}}}_{ij},$ respectively, as functions of concentration and water content as well as membrane equivalent weight and the modulus of the hydrophobic matrix of a dry membrane;^32,33 the required system parameters for these calculations are discussed therein. The properties specific to VRFBs are the infinite dilution diffusion coefficients in water and viscous volume of vanadium ions and vanadium-bisulfate ion pairs in solution, ${{\mathfrak{D}}}_{i,0}^{\infty }$ and ${\tilde{V}}_{i},$ respectively, the vanadium-sulfate binary diffusion coefficient ${{\mathfrak{D}}}_{i,{\mathrm{SO}}_{4}}$, the vanadium-sulfate and vanadium-membrane interaction parameters ${\beta }_{i,{{\rm{SO}}}_{4}}$ and ${\beta }_{i,{\rm{M}}}$, ${\rm{V}}({\rm{III}})$- and ${\rm{V}}({\rm{IV}})$-bisulfate association constants ${K}_{i-{\mathrm{HSO}}_{4}}$, and solvent/vanadium binding constant ${k}_{i}.$

Unfortunately, there are relatively few measurements of vanadium thermodynamic and transport properties in all its oxidation states at well-defined conditions. As described in Supplemental Material (SM), reported experimental values are used when available, and otherwise, parameters of cations of similar charge number approximate those of the vanadium ions. For ${\rm{V}}({\rm{IV}})$ and ${\rm{V}}({\rm{V}})$ we fit ${\beta }_{i,{\rm{M}}}$ to measured vanadium, sulfate, and water uptake in the membrane as a function of sulfuric acid and vanadium sulfate concentration detailed in Part I.³² Similarly, for ${\rm{V}}({\rm{III}}),$ we fit ${\beta }_{i,{\rm{M}}},$ ${\beta }_{i,{{\rm{SO}}}_{4}}$ and ${k}_{i}$ to isotherm measurements. These parameters are plausible values for this system; however, they are estimates. Consequently, there is a strong need for fundamental thermodynamic studies of vanadium-ion properties in their different oxidation states (i.e. activity and osmotic coefficients as a function of concentration).

Throughout this paper, we will consider a reference VRFB system containing electrodes at a fixed composition and a membrane with fixed properties, which are provided in Table II. Calculations are made for this system unless stated. The electrolyte concentrations match those in experimental cells as described in Refs. 1, 3, 51 unless otherwise stated. In those studies, an initial solution of 1.5 mol dm⁻³ vanadium IV sulfate in 2.6 mol dm⁻³ sulfuric acid was twice charged with the positive electrolyte replaced after each charge. As outlined in SM, after this charging process and at a state of charge of 50%, the negative electrolyte has 0.75 mol dm⁻³ of both V(II) and V(III) and a total sulfate concentration of 4.66 mol dm⁻³. The positive electrolyte has the same concentrations of V(IV) and V(V) and a total sulfate concentration of 4.11 mol dm⁻³. The electrolyte solution mass density was 1.19 g cm⁻³.⁵¹ The separator in the reference system is a Nafion 212 membrane with a dry-membrane thickness, ${l}_{{\rm{M}}}^{0}$, of 51 μm. The membrane equivalent weight, $EW$ (mass-polymer per mole-sulfonate group) is 1100.

Table II.  Properties of the reference VRFB system.

Properties	Parameter	Value
Negative electrode composition	${c}_{{\rm{V}}\left({\rm{II}}\right)}^{\beta }$	0.75 mol dm−3
	${c}_{{\rm{V}}\left({\rm{III}}\right)}^{\beta }$	0.75 mol dm−3
	${c}_{{{\rm{SO}}}_{4}}^{\beta }$	4.66 mol dm−3
Positive electrode composition	${c}_{{\rm{V}}\left({\rm{IV}}\right)}^{{\epsilon }}$	0.75 mol dm−3
	${c}_{{\rm{V}}\left({\rm{V}}\right)}^{{\epsilon }}$	0.75 mol dm−3
	${c}_{{{\rm{SO}}}_{4}}^{{\epsilon }}$	4.11 mol dm−3
Polymer backbone modulus	${E}_{b}^{0}$	154 MPa
Membrane equivalent weight	$EW$	1100 g-polymer/mol SO3

Other membrane properties considered are hydrophobic matrix modulus of the dry membrane, ${E}_{{\rm{b}}}^{0},$ Archie's tortuosity scaling parameter, $\chi ,$ hydrophilic domain geometric transport factor, $G,$ and spacing between hydrophilic domain in dry membrane, ${d}^{0}$. $EW$ sets the intrinsic concentration of ions in the membrane absent co-ions, ${E}_{{\rm{b}}}^{0}$ limits the extent the membrane can swell with water from the surrounding environment, $G$ affects friction between the membrane and absorbed water and aqueous ions, and ${d}^{0}$ and $\chi$ dictate the hydrophilic domain size and network tortuosity, respectively, at a given water volume fraction. Unless otherwise stated, membrane properties are those of the Nafion PFSA separator and are detailed in Parts I and II.^32,33

Figure S1 (available online at stacks.iop.org/JES/167/013549/mmedia) outlines calculation of membrane properties, fluxes, and performance metrics. For a specified electrolyte composition in the negative and positive electrodes, Eqs. 14 and 15 specify ${x}_{i}$ at $\delta ^{\prime}$ and $\delta ^{\prime\prime}$. Part I outlines calculation of ${\mu }_{i}^{{\rm{ex}}}.$³² Given the composition of the species in the membrane, Part II outlines calculation of ${{\mathfrak{D}}}_{ij}.$³³ Equation 24 gives ${\chi }_{ij}$ where the Python package Numdifftools V.0.9.39 calculates the Jacobian, $\partial {\mu }_{i,{\rm{M}}}^{{\rm{ex}}}/\partial RT\,\mathrm{ln}({x}_{j}/{x}_{{\rm{M}}}^{{z}_{j}/{z}_{{\rm{M}}}})$ using a first-order, forward-finite-difference approximation with a step size of $(1-{x}_{0})/100$ for water and ${x}_{i}/100$ for ions. Given ${\chi }_{ij}$ and ${{\mathfrak{D}}}_{ij},$ Eq. 21 calculates ${\boldsymbol{D}}.$ The total concentration ${c}_{{\rm{T}}}$ is set to 55.3 mol dm⁻³. For ideal volumetric mixing and isotropic swelling, the membrane thickness is ${l}_{{\rm{M}}}={l}_{{\rm{M}}}^{0}{(1+\lambda {\bar{V}}_{0}{\rho }_{{\rm{M}}}^{0}/EW)}^{1/3},$ $\lambda$ is the average water content ($={x}_{0}/{x}_{{\rm{M}}}$), ${\bar{V}}_{0}$ is the partial molar volume of water (=18.1 cm³ mol⁻¹), and ${\rho }_{{\rm{M}}}^{0}$ is the mass density of the dry membrane (=2.1 g cm⁻³).

To solve Eq. 19 for each species, we extend to electrolyte systems a method outlined by Krishna et al.^52,53 for uncharged systems. Using this approach, Appendix B shows that for a specified current density, the species fluxes are specified by iterating over a set of six transcendental algebraic equations and neglecting changes across the membrane of: the thermodynamic factor ${{\rm{\chi }}}_{ij}$ because the solution composition from the most abundant species (sulfate, the membrane, water, and protons) changes relatively little across the membrane; ${c}_{{\rm{T}}}$ because the molar density of the membrane change only slightly with composition; ${{\mathfrak{D}}}_{ij}$ because the dominant factors its depends upon--ionic strength and water content inside the membrane--are relatively constant; and $\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}/\partial z$ because the membrane concentration changes little between the electrodes. We evaluate ${{\mathfrak{D}}}_{ij}$, ${\chi }_{ij}$, and λ using the mean of the composition of the membrane at the interface with the two electrodes ${x}_{i}^{\bar{\delta }}$ ($=({x}_{i}^{{\delta }^{\prime }}+{x}_{i}^{{\delta }^{\prime\prime }})/2$), calculating the fraction of vanadium and protons that form ion pairs with sulfate at this composition. SM shows that the error introduced by these assumption is small. The derivative of $\mathrm{ln}\,{x}_{{\rm{M}}}$ is approximated as $\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}/\partial z\approx \,\mathrm{ln}({x}_{{\rm{M}}}(z)/{x}_{{\rm{M}}}(z=0))/{l}_{{\rm{M}}}.$ In Results and Discussion, we show this assumption is also reasonable.

Upon specifying species fluxes at a given charge and discharge current density, Eq. 11 gives ${i}_{{\rm{x}}}$ and, using Eq. 12, ${\varepsilon }_{{\rm{q}}}.$ Integration of Eq. 27 gives ${\rm{\Phi }}.$ ${\rm{\Phi }}$ is a linear function of $z$ for the assumption of constant ${{\mathfrak{D}}}_{ij}.$ Upon calculating ${\rm{\Delta }}{\rm{\Phi }},$ Eqs. 6, 2, and 1 give $V,$ ${\varepsilon }_{{\rm{v}}}$, ${\varepsilon }_{{\rm{e}}}$, and ${{\rm{\Psi }}}_{{\rm{d}}}.$ To find the maximum discharge cell power density ${{\rm{\Psi }}}_{{\rm{d}}}$ at ${\varepsilon }_{{\rm{e}}}=80 \%$ (our adopted design criteria), we iteratively solve these equations for varying $i$ until finding the maximum.⁵⁴ As previously stated, the electrode composition and ${N}_{i}$ are assumed constant during charge or discharge. This condition is achieved for incremental changes in $SOC$ during cycling.

We performed sample-based sensitivity analysis for the impact of membrane properties $EW,$ ${E}_{{\rm{b}}}^{0}$, $\chi ,$ $G,$ and ${d}^{0}$ on transport properties and cell performance. SM details the sampling technique and range for these properties. The sensitivity analysis uses a cell current density of 200 mA cm⁻² and the electrolyte composition is that of the reference conditions (see Table II).

Representative VRFB Membrane Transport

Transport properties

Equation 21 shows that the diffusion matrix ${\boldsymbol{D}}$ captures transport coupling between all species due to thermodynamic interactions that ${\boldsymbol{\chi }}$ specifies and for and frictional interactions that ${{\boldsymbol{B}}}^{-1}$ incorporates. For a thermodynamically ideal solution, ${\chi }_{ii}=1$ and ${\chi }_{ij}=0.$ In a dilute solution (${x}_{i}\to 0$ and ${x}_{{\rm{M}}}\to 1$), the matrix ${\boldsymbol{D}}$ has the species/membrane diffusion coefficients on the diagonal while off-diagonal elements are equal to zero.

For Nafion in a VRFB, Fig. 2a plots the mean calculated diffusion coefficients ${D}_{ij}$ on a heat map colored from dark blue (large positive coefficients) through white (coefficient equal to zero) to dark red (large negative coefficient) and the thermodynamic factor ${\chi }_{ij}$ is given in parenthesis. The calculations show that transport in this system is far more complex than diffusion in an ideal, dilute solution. Although some off-diagonal terms of the diffusion matrix are small, they are not categorically negligible and in some cases are larger than the on-diagonal terms. For example, a concentration gradient of ${\rm{V}}\left({\rm{IV}}\right){{\rm{M}}}_{2}$ causes more than twice the flux of ${\rm{V}}({\rm{V}})$ compared to an equal concentration gradient in ${\rm{V}}\left({\rm{V}}\right){\rm{M}}.$ A few of the diffusion coefficients are negative, indicating that the flux of one species will take place up the concentration gradient of another, holding all else constant. In particular, for all vanadium species $i,$ ${D}_{{{\rm{SO}}}_{4}i}$ is negative and large in magnitude.

Zoom In Zoom Out Reset image size

Thermodynamic nonidealities contribute to transport coupling in particular. ${\chi }_{ij}$ for vanadium species $i$ and $j$ are close to ideal because positively charged ions do not significantly change the excess chemical potential of other positively charged species. ${\chi }_{ij}$ for vanadium/water and /sulfate are mostly large because of the significant effect interactions between these species have on the excess free energy of the system.

Figure 2b shows a heat map of charge number-normalized transference numbers. Absent concentration gradients, protons will carry most of the current because ${t}_{{\rm{H}}}=0.94.$ Sulfate ions are the second largest carrier, and vanadium species each transport about 3% of the current.

Transport at open circuit

The concentration difference between electrodes creates concentration gradients in the membrane. Figure 3 shows the mole-fraction profile of each species in the membrane between the interfaces with the negative electrode ($z=0$) and the positive electrode ($z={l}_{{\rm{M}}}$) calculated using the procedure outlined in Appendix B with the no current density (i.e. $i=0$). For transport that is not coupled (i.e. dilute-solution approximation), the mole-fractions of the species will decrease linearly between the two electrodes. This behavior is not present in this system, indicating that transport coupling plays an important role.

Zoom In Zoom Out Reset image size

The vanadium species mole fractions in the negative electrode (V(II) and V(III)) goes to zero on approaching the positive electrode-interface, and the vanadium species mole fractions in the positive electrode (V(IV) and V(V)) decrease towards the negative electrode/membrane interface. The mole fractions of protons, sulfate, and membrane fixed-charged groups are relatively constant across the membrane with the latter two showing a slight decrease/increase towards the negative electrode and having a maximum/minimum away from the electrode interfaces. The dotted line in Fig. 3 shows $\mathrm{ln}({x}_{{\rm{M}}}(z)/{x}_{{\rm{M}}}(z=0))/{l}_{{\rm{M}}}$, which generally matches $\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}/\partial z$ and illustrates that the assumption equating the two is reasonable.

The water mole fraction decreases from the negative electrode to the positive electrode/membrane interface. Although the water mole fraction of the two electrode electrolytes are nearly the same, ${{\rm{\Gamma }}}_{0}$ at the negative electrode is larger than at the positive electrode causing a concentration gradient inside the membrane.

The concentration gradients (see Fig. 3) drive fluxes according to the matrix of diffusion coefficients (see Fig. 2). Using ${\rm{V}}({\rm{II}})$ as an example, Fig. 4 plots the contribution to ${N}_{{\rm{V}}\left({\rm{II}}\right)}$ (see Eq. 19) from the concentration gradient of species $j$ referenced to the membrane (i.e. ${D}_{{\rm{V}}({\rm{II}})j}\partial {c}_{j,{\rm{M}}},$ where we abbreviate $\partial /\partial z$ as $\partial$) as a function of position across the membrane. Each species plays an important role in driving transport with their contribution varying across the membrane. Where a concentration gradient is large, that species drives more transport. For example, around $z=45$ μm, vanadium gradients are largest and contribute substantially to ${\rm{V}}\left({\rm{II}}\right)$ transport. Approaching the positive electrode interface ($z=$ 58 μm), ${x}_{{\rm{V}}\left({\rm{II}}\right)}\to 0,$ and leading to ${D}_{{\rm{V}}\left({\rm{II}}\right)j\ne {\rm{V}}({\rm{II}})}\to 0;$ as this point, ${\rm{V}}({\rm{II}})$ transport is then solely caused by the ${\rm{V}}({\rm{II}})$ concentration gradient.

Zoom In Zoom Out Reset image size

The water concentration gradient is a key transport driving force. Although the water-driven diffusion coefficient ${D}_{{\rm{V}}({\rm{II}})0}$ is not large (see Fig. 2, Row 3, Column 1), the water concentration gradients are (see Fig. 3), thereby causing high water-related fluxes. The term $\partial {c}_{{{\rm{SO}}}_{4}{{\rm{M}}}_{-2}}$ ($=\partial {c}_{{\rm{SO}}4}-\displaystyle \frac{{c}_{{{\rm{SO}}}_{4}}{z}_{{\rm{SO}}4}}{{c}_{{\rm{M}}}{z}_{{\rm{M}}}}\partial {c}_{{\rm{M}}}$) also contributes significantly to transport. Although $\partial {c}_{{\rm{SO}}4}$ is relatively small, as Fig. 3 shows, $-\displaystyle \frac{{c}_{{{\rm{SO}}}_{4}}{z}_{{\rm{SO}}4}}{{c}_{{\rm{M}}}{z}_{{\rm{M}}}}\partial {c}_{{\rm{M}}}$ is large and positive and ${D}_{{\rm{V}}\left({\rm{II}}\right){{\rm{SO}}}_{4}}$ is negative. As a consequence, the term $-{D}_{{\rm{V}}\left({\rm{II}}\right){{\rm{SO}}}_{4}}\partial {c}_{{{\rm{SO}}}_{4}{{\rm{M}}}_{-2}}$ drives ${\rm{V}}({\rm{II}})$ towards the positive electrode. Moreover, at the negative electrode ($z=0$), $\partial {c}_{{\rm{V}}({\rm{II}})}$ is small but $-\displaystyle \frac{{c}_{{\rm{V}}({\rm{II}})}{z}_{{\rm{V}}({\rm{II}})}}{{c}_{{\rm{M}}}{z}_{{\rm{M}}}}\partial {c}_{{\rm{M}}}$ and ${D}_{{\rm{V}}\left({\rm{II}}\right){\rm{V}}({\rm{II}})}$ are large and negative. The result is that $-{D}_{{\rm{V}}\left({\rm{II}}\right){\rm{V}}({\rm{II}})}\partial {c}_{{\rm{V}}({\rm{II}}){{\rm{M}}}_{2}}$ pulls vanadium towards the negative electrode.

Figure 5 shows the calculated normalized mean flux contribution of each species $i$ due to the concentration gradients of each other species $j$ referenced to ${\rm{M}}$ (i.e. ${\int }^{}{D}_{ij}\partial {c}_{j,{\rm{M}}}dz/\sum _{k\ne {\rm{M}},n}| {\int }^{}{D}_{ik}\partial {c}_{k,{\rm{M}}}dz|$). The left side of Fig. 5 gives the total flux of each species multiplied by $F.$ The points show species fluxes for an thermodynamically ideal solution (i.e. ${{\rm{\chi }}}_{ij}=0$ for $i\ne j$ and ${{\rm{\chi }}}_{ii}=1$).

Zoom In Zoom Out Reset image size

As the large magnitude off-diagonal elements in Fig. 2 show, concentration gradients of species $i$ causes transport of species $j.$ This finding is further evidence that dilute-solution theories in which a species transport is driven by only its own concentration gradient (i.e. ${N}_{i}=f(\partial {c}_{i})$) is not sufficient to calculate transport accurately. To quantify this error, the term in parenthesis on the left side of Fig. 5 shows that the dilute-solution approximation (i.e. that the only friction on a species is due to the membrane, ${{\mathfrak{D}}}_{ij\ne {\rm{M}}}\to \infty$) predicts fluxes that, on average (i.e. mean), deviate from the concentrated-solution model by 772%.

For all species, the principle transport driving force is the water concentration gradient. Water drives positively charged species to the positive electrode, resulting in a curvature of the vanadium species concentration gradients towards the positive electrode. To maintain electroneutrality, the water concentration gradient pushes sulfate to the negative electrode.

Mass-transport of ions and water is coupled because both the frictional interactions between species are nonzero and the excess chemical potential of species depends on the concentration of other species. When describing the solution as thermodynamically ideal (circles in Fig. 5), the behavior of the coupling between species is observed, since the thermodynamic nonidealities lead to quantitative differences from the ideal-solution approximation.

In practice, an operating cell will undergo numerous cycles and electrolyte rebalancing such that the composition in the electrodes varies, changing the concentration gradients. Consequently, the specific direction and contributions to fluxes in Fig. 5 do not apply throughout VRFB operation. Incorporating concentrated-solution theory into a fully transient model of VRFB operation would give insight into how these driving forces change during cell operation.

Transport under applied potentials

Under an applied potential, charged species migrate with the current. Figure 6 shows the flux of each vanadium species $i$ (multiplied by ${z}_{i}F$) as a function of current density $i.$ When current is positive, the concentration gradients acting on V(II) and V(III) are aligned and increasing current linearly increases vanadium flux. At negative currents, migration and diffusion are opposed causing the fluxes of V(II) and V(III) to go to zero. V(V) and V(IV) follow opposite scenarios than V(II) and V(III). This description is qualitatively consistent with previous work using a dilute-solution theory framework. For these electrolyte concentrations, the diffusional contributions causes a net flux of vanadium towards the positive electrode when the current is zero as shown in Fig. 5.

Zoom In Zoom Out Reset image size

Figure 7 shows the mean contributions outlined in Eq. 19 to the V(II) flux ${N}_{{\rm{V}}\left({\rm{II}}\right)}$ from each diffusion term, ${D}_{{\rm{V}}\left({\rm{II}}\right)j}\partial {c}_{j,{\rm{M}}},$ (see Fig. 5) and migration, $i{t}_{{\rm{V}}\left({\rm{II}}\right)}/{z}_{{\rm{V}}\left({\rm{II}}\right)}$ as a function of current density. As the current density increases, the fraction of flux due to migration increases. The contributions to flux from diffusion change slightly as the current density increases because migration alters the concentration gradients across the membrane. Above 250 mA cm⁻², the majority of the vanadium V(II) flux is due to migration.

Zoom In Zoom Out Reset image size

Figures 4, 5, and 7 show that coupling plays an important role in transport though VRFB membranes. The contributions from water concentration gradients are particularly large.⁴⁴ The calculations show that in cases that require only semi-quantitative predictions, neglecting the terms ${D}_{ij}\partial {c}_{j,{\rm{M}}}$ for $j\ne i,0$ will cause ∼25% error to the calculated fluxes at OCP with the error decreasing with increasing current density. As a result, mathematical models and experimental analysis that rely on extended versions of dilute-solution theories that incorporate the effect of water on vanadium fluxes will be qualitatively correct in many cases.

Cell performance

Figure 8 shows that the calculated (lines) and measured³ (symbols) voltage and Coulombic efficiencies, ${\varepsilon }_{{\rm{v}}}$ and ${\varepsilon }_{{\rm{q}}},$ respectively, as a function of current density. The calculated and measured efficiencies correspond to slightly different scenarios. The model calculates performance at a constant 50% $SOC$ (i.e. transport coefficients and concentration gradients are constant), whereas in the experiments charge and discharge between two voltage windows corresponding to different $SOC$s. Despite these slightly different conditions, the model and experiment are in good agreement.

Zoom In Zoom Out Reset image size

As the current density increases, the Coulombic efficiency generally increases as the effects of crossover are reduced in comparison to the energy extracted from the cell. However, rising current density decreases voltage efficiency because the membrane and cell incur larger ohmic losses. As a consequence, the energy density increases initially, reaches a maximum, and then decreases with increasing current density. In the absence of resistance from the membrane ${R}_{{\rm{M}}}=0$ (i.e. ${\rm{\Delta }}{\rm{\Phi }}=0$), ${\varepsilon }_{{\rm{v}}}$ is higher, but non-membrane resistances limit the benefits of ultra-high conductive membranes.

Figure S2 shows voltage and Coulombic efficiency calculated using dilute-solution theory that only accounts for friction with the membrane (i.e. ${{\mathfrak{D}}}_{ij\ne {\rm{M}}}\to \infty$) and an extended dilute-solution theory that also accounts for frictions with water (i.e. ${{\mathfrak{D}}}_{ij\ne 0,{\rm{M}}}\to \infty$). The dilute-solution theory severely miscalculates cell efficiencies. For Coulombic efficiency, the extended-dilute-solution theory closely agrees with concentrated-solution theory but slightly overestimates voltage efficiency. Based on this analysis, researchers should account for the effect of coupled-transport modes and, in particular, the role of ion-water interactions on cell performance.

The molecular transport and thermodynamic models presented in Parts I and II calculate properties as a function of external electrolyte composition and membrane chemistry.^32,33 The following sections use these models to determine how changes away from the reference system (see Table II) to operating conditions and membrane properties affect performance.

Impact of electrolyte concentration

Concentration of species in the external electrolyte solution strongly impact membrane transport properties. Figure 9 shows membrane conductivity $\kappa$ (a), V(II)0 diffusion coefficient ${D}_{{\rm{V}}\left({\rm{II}}\right)0}$ (b), and V(II) transference number ${t}_{{\rm{V}}({\rm{II}})}$ (c) calculated at membrane composition ${x}_{i}^{\bar{\delta }}$ as a function of the total vanadium molarity in the electrodes, and the lines are varying arithmetic average sulfate concentration in the electrodes. Figure S3 shows the membrane water content (a), proton molality (b), and V(II) molality (c) at the same conditions. For brevity, we plot the properties of V(II), which are representative of the transport properties of the other vanadium species.

Zoom In Zoom Out Reset image size

As the concentration of vanadium in the external electrolyte rises, the proton content in the membrane decreases, as Fig. S3b shows. The loss of protons, which have high mobility, decreases membrane conductivity. As Fig. S3c shows, the vanadium content in the membrane also increases with rising vanadium concentration in the electrolyte leading to larger vanadium transport properties ${D}_{{\rm{V}}\left({\rm{II}}\right)0}$ and ${t}_{{\rm{V}}\left({\rm{II}}\right)}.$ Adding sulfuric acid increases the membrane proton concentration. Higher proton concentrations displace vanadium from the membrane, as Figs. S4b and c show. Consequently, at high vanadium concentrations in the electrode, rising sulfate concentration increases conductivity and decreases vanadium transport properties, as Fig. 9 shows. However, adding sulfuric acid in the external electrolyte also dehydrates the membrane, as Fig. S3a shows and Part I discusses.³² At low vanadium concentration in the electrodes, the dehydrating effects of additional sulfate dominates, causing membrane conductivity to decrease.

Figure 10 shows the calculated Coulombic and voltage efficiencies, ${\varepsilon }_{{\rm{q}}}$ and ${\varepsilon }_{{\rm{v}}}$ respectively, (a) and discharge power density ${{\rm{\Psi }}}_{{\rm{d}}}$ (b) at the current density that maximizes ${{\rm{\Psi }}}_{{\rm{d}}}$ at ${\varepsilon }_{{\rm{e}}}$ =80% (termed optimal power density) as a function of the vanadium concentration and for different mean sulfate concentrations in the electrodes. To show how sensitive the optimal power density is to the target energy efficiency, shaded regions in Fig. 10b show power density at ±0.2% of the target energy density. As Fig. 9 shows, rising vanadium concentration in the electrode electrolyte increases the vanadium transference numbers and diffusion coefficients, decreasing Coulombic efficiency. To maintain an 80% energy efficiency, the current density decreases to lower ohmic losses and increase voltage efficiency. The cumulative effect lowers the optimal power density with increasing vanadium concentration. Moreover, the lower Coulombic efficiency at high vanadium concentration makes the optimal power density less sensitive to the required energy efficiency (i.e. narrower shaded region in Fig. 9).

Zoom In Zoom Out Reset image size

Increasing sulfuric acid concentration in the electrolyte (shown in Fig. 10 as moving from dot-dashed to dashed to solid lines) exchanges vanadium in the membrane for protons, increasing membrane conductivity and decreasing vanadium transport coefficients, as Fig. 9 shows. The larger resulting Coulombic efficiency allows the system to incur more ohmic losses by operating at higher current densities while maintaining the target energy efficiency. As a result, increasing sulfuric-acid concentration in the electrolyte increases the optimal power density. However, at low vanadium concentrations, increasing sulfate concentration lowers membrane conductivity, as Fig. 9 shows. The addition of sulfuric acid will therefore be less effective at low active-species concentrations. Moreover, the improved optimal power density at low vanadium and high sulfuric-acid concentrations may not result in higher return on investment for VRFB systems if balance-of-plant costs increase due to lower energy density or higher pumping costs.⁴²

Impact of membrane properties

Figure 11 shows membrane conductivity $\kappa$ (a), V(II)/0 diffusion coefficient ${D}_{{\rm{V}}\left({\rm{II}}\right)0}$ (b), and V(II) transference number ${t}_{{\rm{V}}({\rm{II}})}$ (c) calculated at membrane composition ${x}_{i}^{\bar{\delta }}$ as a function of EW and different dry-membrane moduli ${E}_{{\rm{b}}}^{0}.$ Figure S4 shows the membrane water content $\lambda$ (a), proton molality ${m}_{{\rm{H}}}^{\bar{\delta }}$ (b), and V(II) molality ${m}_{{\rm{V}}({\rm{II}})}^{\bar{\delta }}$ (c) at the same conditions.

Zoom In Zoom Out Reset image size

Increasing EW decreases membrane water content, as Fig. S4a shows, thereby decreasing hydophilic water volume fraction, increasing tortuosity of the hydophilic domains, and decreasing the pore size. As a result, membrane conductivity $\kappa$ and V(II)/water diffusion coefficient ${D}_{{\rm{V}}\left({\rm{II}}\right)0}$ decrease, as Figs. 11a and 11b show. Furthermore, increasing EW decreases the concentration of vanadium and protons in the membrane because there are fewer oppositely charged sulfonate groups for these ions to interact with. As EW increases, the vanadium concentration decreases more rapidly than proton concentration in the membrane, as Fig. S4b and c show, decreasing the vanadium transference number, as Fig. 11c shows.

Figure 11 shows (moving from dot-dashed to solid lines) that increasing the membrane modulus impacts transport properties similary to increasing EW. Raising the membrane modulus decreases water content, as Fig. S4a shows. The lower water content of the higher modulus membranes results in lower membrane conductivity and vanadium-water diffusion coefficients. As with EW, increasing the membrane modulus decreases the vanadium content in the membrane more rapdily than the proton content, decreasing the vanadium transference number.

Figure 12 shows the calculated Coulombic ${\varepsilon }_{{\rm{q}}}$ and voltage ${\varepsilon }_{{\rm{v}}}$ efficiencies (a) and discharge power density ${{\rm{\Psi }}}_{{\rm{d}}}$ (b) at the current density that maximizes ${{\rm{\Psi }}}_{{\rm{d}}}$ at ${\varepsilon }_{{\rm{e}}}$ =80% as a function of equivalent weight EW and dry-membrane modulus ${E}_{{\rm{b}}}^{0}.$ As Fig. 11 shows, increasing EW decreases the vanadium transference numbers and diffusion coefficients, increasing Coulombic efficiency. However, conductivity also falls with increasing EW, lowering voltage efficiency and cell potential. The benefits of higher Coulombic efficiency do not outweigh the costs to voltage efficiency and cell potential. The net effect of increasing EW is a lower optimal power density.

Zoom In Zoom Out Reset image size

Decreasing the membrane modulus (shown in Fig. 12 as moving from dot-dashed, to dashed, to solid lines) increases membrane conductivity, vanadium diffusion coefficients, and vanadium transference numbers, as Fig. 11 indicates. At high EWs, the resulting increase in cell potential and voltage efficiency more than compensates for lower Coulombic efficiency, leading to a higher optimal power density. These findings agree with experimental studies of cell performance with different EW membranes and varying mechanical reinforcement.^6,9,30 At low EWs (<900 g-polymer/mole SO₄), ameliorating the poor Coulombic efficiency by increasing membrane modulus offsets the voltage-efficiency losses and results in highest optimal power density.

The model results show that transport properties in VRFB membranes (e.g. conductivity, diffusivities, etc) result from a complex interplay of multiple phenomena that are mediated by structural characteristics of the membrane (e.g. membrane modulus and EW). These transport properties, in turn, dictate VRFB performance (e.g. Coulombic efficiency). To provide general structure-property-performance relationships, we perform sample-based sensitivity analysis. This analysis (see Numerical Implementation) calculates properties and the performance of membranes over a range of equivalent weights, moduli, tortuosities, and hydrophilic domain shapes and sizes. This method is analogous to synthesizing a thousand membranes with attributes spanning this parameter space and measuring their properties and in-cell performance.

Figures S5, S6, and S7 show pair plots of structural attributes ($EW,$ ${E}_{{\rm{b}}}^{0},$ $\chi ,$ $G,$ and ${d}^{0}$), exemplar transport and uptake properties ($\kappa ,$ ${t}_{{\rm{II}}},$ ${D}_{0{\rm{II}}},$ ${D}_{{\rm{IIII}}},$ ${m}_{{\rm{II}}},$ ${m}_{{\rm{H}}},$ and $\lambda$), and cell performance at a current density of 200 mA cm⁻² (${{\rm{\Psi }}}_{{\rm{d}}},$ ${\varepsilon }_{{\rm{q}}},$ and ${\varepsilon }_{{\rm{v}}}$) and reference conditions (see Table II). Figure S8 gives the least-squares linear regression coefficients between the natural log of each variable. Based on the sample-based sensitivity analysis, Fig. 13 shows a chord diagram of correlation between (a) structural attributes of the membrane and transport and uptake properties and (b) membrane properties and cell performance metrics scaling the chord size by the Spearman rank correlation coefficient ${r}_{{\rm{s}}}$ squared (Fig. S9 provides all values of the Spearman coefficients between variables). Spearman coefficient is a measure of correlated monotonicity between variables.

Zoom In Zoom Out Reset image size

Figure 13a illustrates that dry-membrane domain spacing, ${d}^{0},$ and membrane modulus, ${E}_{{\rm{b}}}^{0},$ are highly correlated to the molality of ions in the membrane and water uptake, respectively, because they directly affect the energetics of membrane partitioning (see Part I).³² The shape and tortuosity of the hydrophilic domains, quantified by $G$ and $\chi ,$ respectively, correlate with transport properties. The domain shape affects the transference number of vanadium, whereas the tortuosity influences the conductivity and diffusion coefficients. ${\rm{EW}}$ correlates with all of these membrane properties, but typically to a lesser extent than the other membrane attributes. Figure 13a suggests which membrane design approaches are likely to influence a property of interest. However, the structure-function correlations are not high (i.e. at most 0.8), indicating that varying the structure of the membrane does not guarantee altered properties. Moreover, structural attributes impact multiple transport properties such that properties cannot be tuned independently.

Figure 13b correlates the membrane properties with in-situ VRFB performance. Membrane conductivity directly governs voltage efficiency, ${\varepsilon }_{{\rm{v}}},$ and discharge power density, ${{\rm{\Psi }}}_{{\rm{d}}},$ at a current density of 200 mA cm⁻² (see Eqs. 1, 6, and 28), giving ${r}_{{\rm{s}}}=1$ between these variables. Because conductivity is highly correlated with vanadium diffusivity, vanadium diffusion coefficients are also correlated with ${\varepsilon }_{{\rm{v}}}$ and ${{\rm{\Psi }}}_{{\rm{d}}}$ but do not have a causal effect (see Fig. S7). The sensitivity analysis shows that Coulombic efficiency is highly correlated with vanadium transference number (${r}_{{\rm{s}}}=-0.98$) and water-vanadium diffusion coefficient (${r}_{{\rm{s}}}=-0.84$) and only weakly correlated with vanadium-vanadium diffusion coefficient (${r}_{{\rm{s}}}=-0.55$). These correlations are consistent with Fig. 7 in showing the crucial roles of migration and water-driven transport in cell operation. Water and ion uptake are weakly correlated to cell performance metrics. As such, screening membranes using ex-situ measurements of conductivity and vanadium transference number are most indicative of cell performance in these conditions.

This work develops a multicomponent, concentrated-solution model of transport in VRFB membranes that accounts for Stefan-Maxwell-Onsager transport couplings between species and thermodynamic nonidealities. The molecular-thermodynamics model outlined in Part I calculates ion and water partitioning into the membrane from the electrolyte in the electrodes and provides thermodynamic factors that influence species diffusion. The molecular-transport model outlined in Part II describes binary-diffusion coefficients between ions, water, and the membrane as a function of water content and composition.^32,33 The resulting multicomponent diffusion coefficient matrix has large off-diagonal elements that dilute-solution theory neglects and that impact the net transport of the species through the membrane. Transport is highly coupled because of both thermodynamic and frictional interactions between species. In particular, water concentration gradients play a dominant role in vanadium crossover. Under an applied current, contributions to net flux from these diffusional modes and migration are aligned or opposed, depending on the direction of transport. Migration dominates at high current densities (>250 mA cm⁻²). The magnitude of coupling suggests that predictions made using dilute-solution theories for VRFB will be misleading. Extended-dilute-solution theories that account for the impact of water chemical potential gradients on vanadium transport are more reliable. Furthermore, membrane-permeation experiments that measure flux under ill-defined concentration gradients will not provide diffusion coefficients that are predictive of transport under other conditions. Even in the case of carefully designed experiments, the magnitude of measured diffusion coefficients should be interpreted in the context of the collective interactions between all species present in the system. Sensitivity analysis indicate that membrane conductivity and vanadium transference number are the best predictors of cell performance metrics.

The cell performance in terms of power density and energy efficiency depends on the complex interplay of the thermodynamic and frictional interactions of all species and their gradients. We studied the effect of sulfate and vanadium concentration in the electrode electrolyte and membrane equivalent weight and modulus on membrane properties, species transport, and cell performance. For a metric of maximum power density at 80% energy efficiency, low-vanadium and high-sulfuric acid concentration electrolytes perform best. However, the solubility of vanadium in sulfuric acid limits the success of this strategy. Membranes with low equivalent weights (i.e. high ion-exchange capacity) perform better. Specifically, decreasing EW from 1500 to 900 improves performance by 20% by increasing membrane conductivity. At high equivalent weights, a low modulus is superior, although perhaps not readily synthesizable.

The work presented here was funded, in part, by the Advanced Research Projects Agency - Energy (ARPA-E), U.S. Department of Energy (DOE) under Award Numbers DEAR0000149 (UTRC) and DEAC0205CH11231 (LBNL with cost share from UTRC). (The information, data, or work presented herein was funded in part by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof). Work by Robert Darling was supported as part of the Joint Center for Energy Storage Research, an Energy Innovation Hub funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under contract DE-AC02–06CH11357. The authors acknowledge fruitful discussions with Michael Gerhardt, Anamika Chowdhury, Ahmet Kusoglu, and Joseph Arthur.

Adding Eq. 18 to Eq. 17 puts the driving force in terms of ionic potential-independent term ${\mu }_{i,{\rm{M}}}$

$$\begin{eqnarray}&&\begin{array}{l}{c}_{i}\displaystyle \frac{\partial {\mu }_{i}}{\partial z}-\displaystyle \frac{{c}_{i}{z}_{i}}{{z}_{{\rm{M}}}}\displaystyle \frac{\partial {\mu }_{{\rm{M}}}}{\partial z}={c}_{i}\displaystyle \frac{\partial {\mu }_{i,{\rm{M}}}}{\partial z}\\ \,=\,\displaystyle \displaystyle \sum _{j\ne i,{\rm{M}}}\displaystyle \frac{RT}{{{ \mathcal D }}_{ij}{c}_{{\rm{T}}}}\left({c}_{i}{N}_{j}-{c}_{j}{N}_{i}\right)-\displaystyle \frac{{c}_{{\rm{M}}}RT}{{c}_{{\rm{T}}}{{ \mathcal D }}_{i{\rm{M}}}}{N}_{i}\\ \,\,\,-\,\displaystyle \frac{{c}_{i}{z}_{i}}{{c}_{{\rm{T}}}{z}_{{\rm{M}}}}\displaystyle \displaystyle \sum _{j\ne {\rm{M}}}\displaystyle \frac{RT}{{{ \mathcal D }}_{j{\rm{M}}}}{N}_{j}\end{array}\end{eqnarray} \tag{ A1 }$$

where ${D}_{iM}={c}_{i}{c}_{j}RT/{c}_{{\rm{T}}}{K}_{ij}.$ For an isothermal, isobaric system, the gradient of ${\mu }_{i,{\rm{M}}}$ is a function of $N-2$ species compositions because the Gibbs-Duhem equation and electroneutrality reduces the degrees of freedom by 2. We expand left side of Eq. A1 using the composition variable $\mathrm{ln}({x}_{j}/{x}_{{\rm{M}}}^{{z}_{j}/{z}_{{\rm{M}}}})$

$$\begin{eqnarray}&&\begin{array}{l}{c}_{i}\displaystyle \frac{\partial {\mu }_{i,{\rm{M}}}}{\partial z}={c}_{{\rm{T}}}{x}_{i}\displaystyle \sum _{j\ne {\rm{M}},n}^{N}\displaystyle \frac{\partial {\mu }_{i,{\rm{M}}}}{\partial \,\mathrm{ln}\left(\tfrac{{x}_{j}}{{x}_{{\rm{M}}}^{{z}_{j}/{z}_{{\rm{M}}}}}\right)}\displaystyle \frac{\partial \,\mathrm{ln}\left(\tfrac{{x}_{j}}{{x}_{{\rm{M}}}^{{z}_{j}/{z}_{{\rm{M}}}}}\right)}{\partial z}\\ \,\,\,\,\,=\,RT{c}_{{\rm{T}}}\displaystyle \sum _{j\ne {\rm{M}},n}^{N}\left[{{\rm{\chi }}}_{ij}\left(\displaystyle \frac{\partial {x}_{j}}{\partial z}-\displaystyle \frac{{z}_{j}{x}_{j}}{{z}_{{\rm{M}}}{x}_{{\rm{M}}}}\displaystyle \frac{\partial {x}_{{\rm{M}}}}{\partial z}\right)\right]\end{array}\end{eqnarray} \tag{ A2 }$$

Substitution of Eq. A2 into Eq. A1 gives

$$\begin{eqnarray}&&\begin{array}{l}{c}_{{\rm{T}}}\displaystyle \sum _{j\ne {\rm{M}},n}^{N}{{\rm{\chi }}}_{ij}\displaystyle \frac{\partial {x}_{j}}{\partial z}=\displaystyle \displaystyle \sum _{j\ne i,{\rm{M}}}\displaystyle \frac{({x}_{i}{N}_{j}-{x}_{j}{N}_{i})}{{{ \mathcal D }}_{ij}}-\displaystyle \frac{{x}_{{\rm{M}}}}{{{ \mathcal D }}_{i{\rm{M}}}}{N}_{i}\\ \,\,-\,\displaystyle \frac{{x}_{i}{z}_{i}}{{z}_{{\rm{M}}}}\displaystyle \displaystyle \sum _{j\ne {\rm{M}}}\displaystyle \frac{{N}_{j}}{{{ \mathcal D }}_{j{\rm{M}}}}+{c}_{{\rm{T}}}\displaystyle \sum _{j\ne {\rm{M}},n}^{N}{{\rm{\chi }}}_{ij}\displaystyle \frac{{z}_{j}{x}_{j}}{{z}_{{\rm{M}}}}\displaystyle \frac{\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}}{\partial z}\end{array}\end{eqnarray} \tag{ A3 }$$

By choosing species $n$ such that it has a charge number different from ${\rm{M}},$ we can eliminate the mole fraction of species $n$ and ${\rm{M}}$ on the right side by using the dependence of mole fractions $\left(1=\sum _{i}{x}_{i}\right)$ and electroneutrality ($\sum _{i}{x}_{i}{z}_{i}=0$), specifically giving ${x}_{n}=({z}_{n}+\sum _{j\ne {\rm{M}},n}{x}_{j}\left({z}_{j}-{z}_{n}\right))/({z}_{n}-{z}_{{\rm{M}}})$ and ${x}_{{\rm{M}}}=({z}_{{\rm{M}}}+\sum _{j\ne {\rm{M}},n}{x}_{j}\left({z}_{j}-{z}_{{\rm{M}}}\right))/({z}_{{\rm{M}}}-{z}_{n}).$ Rearranging A3 into a matrix form gives

$$\begin{eqnarray}&&{\boldsymbol{\chi }}{\boldsymbol{x}}^{\prime} ={\boldsymbol{\Psi }}{\boldsymbol{x}}+{\boldsymbol{\zeta }}\end{eqnarray} \tag{ A4 }$$

where ${\boldsymbol{x}}^{\prime}$ and ${\boldsymbol{x}}$ are $N-2$ vectors with entries $\partial {x}_{i}/\partial z$ and ${x}_{i},$ respectively; ${\boldsymbol{\Psi }}$ is a $N-2$ by $N-2\,$matrix containing entries

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Psi }}}_{ii}=-\displaystyle \frac{{N}_{i}}{{c}_{{\rm{T}}}}\left(\displaystyle \frac{{z}_{i}-{z}_{n}}{{{ \mathcal D }}_{i{\rm{M}}}\left({z}_{n}-{z}_{{\rm{M}}}\right)}+\displaystyle \frac{{z}_{i}-{z}_{{\rm{M}}}}{{{ \mathcal D }}_{in}\left({z}_{{\rm{M}}}-{z}_{n}\right)}\right)\\ \,\,\,\,+\,\displaystyle \displaystyle \sum _{k\ne {\rm{M}},i}\displaystyle \frac{{N}_{k}}{{c}_{{\rm{T}}}{{ \mathcal D }}_{ik}}+\displaystyle \frac{{z}_{i}}{{z}_{{\rm{M}}}}\left({{\rm{\chi }}}_{ii}\displaystyle \frac{\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}}{\partial z}-\displaystyle \displaystyle \sum _{k\ne {\rm{M}}}\displaystyle \frac{{N}_{k}}{{c}_{{\rm{T}}}{{ \mathcal D }}_{k{\rm{M}}}}\right)\end{array}\end{eqnarray} \tag{ A5 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Psi }}}_{ij}={{\rm{\chi }}}_{ij}\displaystyle \frac{{z}_{j}}{{z}_{{\rm{M}}}}\displaystyle \frac{\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}}{\partial z}-\displaystyle \frac{{N}_{i}}{{c}_{{\rm{T}}}}\left(\displaystyle \frac{1}{{{ \mathcal D }}_{ij}}+\displaystyle \frac{{z}_{j}-{z}_{n}}{{{ \mathcal D }}_{i{\rm{M}}}\left({z}_{n}-{z}_{{\rm{M}}}\right)}\right.\\ \,\,\,+\left.\,\displaystyle \frac{{z}_{j}-{z}_{{\rm{M}}}}{{{ \mathcal D }}_{in}\left({z}_{{\rm{M}}}-{z}_{n}\right)}\right)\end{array}\end{eqnarray} \tag{ A6 }$$

for $i\ne j;$ and ${\boldsymbol{\zeta }}$ is a $N-2$ vector containing entries

$$\begin{eqnarray}&&{\zeta }_{i}=-\displaystyle \frac{{N}_{i}}{{c}_{{\rm{T}}}}\left(\displaystyle \frac{{z}_{n}}{{{ \mathcal D }}_{i{\rm{M}}}\left({z}_{n}-{z}_{{\rm{M}}}\right)}+\displaystyle \frac{{z}_{{\rm{M}}}}{{{ \mathcal D }}_{in}\left({z}_{{\rm{M}}}-{z}_{n}\right)}\right)\end{eqnarray} \tag{ A7 }$$

Equations A5–A7 are for $i$ and $j$ not equal to species ${\rm{M}}$.

Krishna et al.^52,53 showed that with boundary conditions on ${x}_{i}$ at $z=0$ and $={l}_{{\rm{M}}}$ the solution to Eq. A4 for constant ${{\rm{\chi }}}_{ij},$ ${c}_{{\rm{T}}},$ and $\tfrac{\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}}{\partial z}$ (we approximate $\tfrac{\partial \,\mathrm{ln}\,{x}_{{\rm{M}}}}{\partial z}$ as $\mathrm{ln}\left.\left(\tfrac{{x}_{{\rm{M}}}\left(z={l}_{{\rm{M}}}\right)}{{x}_{{\rm{M}}}\left(z=0\right)}\right)/{l}_{{\rm{M}}}\right)$ is

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{x}}-{\boldsymbol{x}}\left(z=0\right)=\left[\exp \left({{\boldsymbol{\Theta }}}^{-1}{\boldsymbol{\Psi }}z\right)-{\bf{I}}\right]\\ \,\,{\left[\exp \left({{\boldsymbol{\Theta }}}^{-1}{\boldsymbol{\Psi }}{l}_{{\rm{M}}}\right)-{\bf{I}}\right]}^{-1}\,\left[{\boldsymbol{x}}\left(z={l}_{{\rm{M}}}\right)-{\boldsymbol{x}}\left(z=0\right)\right]\end{array}\end{eqnarray} \tag{ B1 }$$

where ${\bf{I}}$ is the identity matrix, the superscript $-1$ denotes the matrix inverse, and the $\exp$ is the matrix exponential. The concentration gradient vector is then

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{x}}\mbox{'}=\left[{{\boldsymbol{\Theta }}}^{-1}{\boldsymbol{\Psi }}\exp \left({{\boldsymbol{\Theta }}}^{-1}{\boldsymbol{\Psi }}z\right)\right]{\left[\exp \left({{\boldsymbol{\Theta }}}^{-1}{\boldsymbol{\Psi }}{l}_{{\rm{M}}}\right)-{\bf{I}}\right]}^{-1}\,\\ \,\,\,\left[{\boldsymbol{x}}\left(z={l}_{{\rm{M}}}\right)-{\boldsymbol{x}}\left(z=0\right)\right]\end{array}\end{eqnarray} \tag{ B2 }$$

We calculate the matrix exponential by diagonalizing the matrix argument and computing the exponential element-wise of the matrix eigenvalues. Eqs. B1 and B2 give the species concentration and concentration gradients as a function of species fluxes.

Equation B2 In conjunction with Eq. 19 calculate flux. The equations are an implicit set of algebraic equations that calculate flux. Krishna et al.^52,53 showed that the following numerical method specifies flux: (1) calculate ${\boldsymbol{\Psi }}$ for ${\boldsymbol{N}}=0;$ (2) calculate ${\boldsymbol{x}}\mbox{'}$ from Eq. B2 at $z=0;$ (3) use calculated ${\boldsymbol{x}}\mbox{'}$ at $z=0$ to calculate ${\boldsymbol{N}}$ using Eq. 19; (4) recalculate ${\boldsymbol{\Psi }}$ using the calculated ${\boldsymbol{N}};$ (5) repeat steps 2–4 until achieving convergence (here, set to 1 × 10^–3 relative change in flux over an iterations).