A Unit Cell Model of a Regenerative Hydrogen-Vanadium Fuel Cell

In this study, a time dependent model for a regenerative hydrogen-vanadium fuel cell is introduced. This lumped isothermal model is based on mass conservation and electrochemical kinetics, and it simulates the cell working potential considering the major ohmic resistances, a complete Butler–Volmer kinetics for the cathode overpotential and a Tafel–Volmer kinetics near mass-transport free conditions for the anode overpotential. Comparison of model simulations against experimental data was performed by using a 25cm 2 lab scale prototype operated in galvanostatic mode at different current density values (50 − 600Am − 2 ). A complete Nernst equation derived from thermodynamic principles was ﬁtted to open circuit potential data, enabling a global activity coefﬁcient to be estimated. The model prediction of the cell potential of one single charge-discharge cycle at a current density of 400Am − 2 was used to calibrate the model and a model validation was carried out against six additional data sets, which showed a reasonably good agreement between the model simulation of the cell potential and the experimental data with a Root Mean Square Error (RMSE) in the range of 0.3–6.1% and 1.3–8.8% for charge and discharge, respectively. The results for the evolution of species concentrations in the cathode and anode are presented for one data set. The proposed model permits study of the key factors that limit the performance of the system and is capable of converging to a meaningful solution relatively fast (s–min).

Redox flow batteries are considered to be an exceptional candidate for grid-scale energy storage.2][3][4] All-Vanadium Redox Flow Batteries (VRFBs) have been considered a promising system due to the limited impact of cross-contamination.However, they have faced challenges related to cost, scale-up and optimization.Current research is also focused on improvement of electrolyte stability for use over a wider temperature window and concentrations, development of electrode materials resistant to overcharge, and mitigation of membrane degradation. 1,2Cost dependency with regarding to vanadium can be mitigated through utilization of new systems that employ only half of the vanadium. 1Recently, a Regenerative Hydrogen-Vanadium Fuel Cell (RHVFC) based on an aqueous vanadium electrolyte V(V) and V(IV) and hydrogen has been introduced 5 and is illustrated schematically in Figure 1.This system contains a porous carbon layer for the positive electrode reaction, membrane and catalyzed porous carbon layer for the negative electrode reaction.Hydrogen evolution, which is an adverse reaction in VRFBs, is here the main anodic process.During discharge, V(V) is reduced to V(IV) and H 2 is oxidized, while the reverse process occurs during charge and H 2 is stored.The vanadium reaction takes place in the positive electrode (cathode), while the hydrogen reaction occurs in the catalyst layer (CL) of the negative electrode (anode).The redox reactions that occur at the electrodes are presented in Equations 2 and 1, and the overall cell reaction in Equation 3, where the charged species VO 2+ and VO + 2 represent the V(IV) and V(V) oxidation states, respectively; and E • is standard potential with the subscripts ca and an referring to positive and negative side respectively.
6][7] An in depth analysis of the system by means of experimental studies has been reported previously by our group. 5,7Yufit et al. 5 studied the performance of a RHVFC and observed that higher current densities or lower catholyte flow rates produce lower coulombic efficiencies, which can be explained by the mass transport limitations.Hewa Dewage et al. 7 studied the loss mechanisms of a RHVFC, observing crossover of vanadium, which could be adsorbed onto the anode Pt catalyst, and a negligible influence of the hydrogen flow rate.A better performance could be obtained if improved component materials and operating conditions were used.Modelling and simulation is an indispensable tool, saving time and reducing cost.To the best of the authors' knowledge, only one RHVFC model has been developed, as part of a PhD project, 8 which presented a one-dimensional steady-state approach.Common VRFB assumptions were considered for the cathode, [8][9][10] while the Membrane Electrolyte Assembly (MEA) was modeled using the assumptions developed by Bernardi and Verburgge. 11,12Only protons were allowed to cross the membrane, which was assumed (3) (4) (5) (6)   Figure 1.Schematic of the RHVFC and model domains.
) unless CC License in place (see abstract).ecsdl.org/site/terms_useaddress.Redistribution subject to ECS terms of use (see 155.198.12.189Downloaded on 2018-02-26 to IP to be fully hydrated since the catholyte is circulating through the cathode that is in direct contact with the membrane. 8A similar assumption was made in a H 2 /Br 2 battery model recently introduced by You et al., 13 but the fully hydrated assumption was only considered at the membrane-cathode interface. 5][16][17][18] Membrane water transport in PEMFC has been addressed following two mechanisms: 14,16,19 diffusive 20 and hydraulic. 11,12The diffusive approach has been the preferred one since it describes more accurately non-flooded membranes. 16Phase change processes such as membrane adsorption/desorption [21][22][23][24] and condensation/evaporation 22,[24][25][26][27] have been also considered in PEMFC models. 16,19Condensation and evaporation have been formulated using a switch function, as well as with overall rate constants.These constants were selected by He et al. 25 to practically generate an instantaneous process while other works have used values from 1 to 10 4 s −1 . 16,19Adsorption and desorption have been formulated considering the departure from the equilibrium water content and to be proportional to rate constants that depend on the water content and temperature. 16,199][30][31][32] It has been observed in VRFBs that vanadium moves across the membrane due to the concentration gradient, while migration enhances or decreases the crossover flux depending on the relative direction of the concentration and potential gradients. 33The crossover flux of VO 2+ has been seen to approach the diffusive flux at Open Circuit Potential (OCP) under large negative current densities, whereas it can reach the sum of migration and electro-osmotic fluxes at large positive current densities. 33A model that explained the crossover of all ionic species in VRFBs has been recently reported, 30,31 which indicated that convection transport produces the major contribution to vanadium ions crossover.For a Br 2 /H 2 cell, minimal crossover during OCP has been observed, along with a strong dependency of crossover of water and bromine species with current density and minimal crossover of hydrogen which does not contribute importantly to the self-discharge of the cell. 33,34In general, it has been stated for redox flow batteries that lower crossover of water and ionic species occurs for thicker membranes and an increase in crossover takes place with higher current densities. 33,34This indicates that the relative importance of migration and electro-osmotic fluxes over diffusion flux increases as both the membrane thickness and current density increase.
Electrochemical reaction kinetics in PEMFCs are not fully understood since the oxygen and hydrogen reactions can be represented by different reaction paths with more than one step. 16,35,36A Butler-Volmer relation has been widely used to model the kinetics, which uses transfer coefficients with values of about 0.5 for both electrodes and neglects the dependency on proton concentration and activity. 14,16This assumption corresponds to the electro-neutrality condition assumed in the polymer electrolyte, 14,16 providing that the proton concentration is equal to the concentration of fixed charge sites.6][37] The Open Circuit Potential (OCP) has been considered to be dependent on the pressure of hydrogen, oxygen and liquid water. 14[41] This equation assumes that the activity coefficient (γ) for all species is equivalent to one, and includes the effect of the cathode proton concentration and Donnan potential, which represents the potential jump at the electrolyte-membrane interfaces due to different proton concentrations in the anolyte and catholyte. 38However, the inconsistency with thermodynamics of this equation has been stated by Pavelka et al., 42 who have introduced a OCP equation for VRFBs derived based on thermodynamic principles (Equation 5).
3][44] The relation proposed by Pavelka et al. underestimated the OCP data when unity activity coefficients are considered. 42o improve the agreement non-ideal activity coefficients should be considered, which can be evaluated from experimental OCP data.The reaction kinetics of VRFBs is commonly modeled using the Butler-Volmer equation, 18 which neglects the effect of the concentration of protons in the positive side. 10,30,41,45It has been recognized that when the transfer coefficients in the Butler-Volmer equation are equal to 0.5, the ability to fit polarization data can be limited. 30,41he incorporation of detailed models in monitoring tools or stack models is currently not feasible. 46Lumped parameter models are suitable for these applications since they are typically based on a small set of differential algebraic equations, which can be fairly simple to integrate and can be accurate enough to describe the battery dynamics. 18umped models consider the time dependent electrochemical phenomena, but assuming a uniform spatial distribution of species concentrations.A VRFB lumped model was developed by Shah et al., 46 which considered mass conservation and the electrochemical processes, and could capture the relation between performance and key properties. 18Additional effects were later included such as thermal processes, 47,48 self discharge, 29,49,50 pumping power and shunt current losses. 51,52A PEMFC lumped model was developed by Pukrushpan et al., 53 which allowed the evolution of component partial pressures and membrane humidity to be studied. 54Effects such as reactant starvation and flooding were later considered to describe the effect of water accumulation on the cell behavior. 55,56n this paper, a simplified model for a unit cell of a hydrogen vanadium system is presented, which was based on the physicochemical phenomena explained by mass conservation, transport mechanisms and electrochemical processes, but maintaining simplicity to allow its use in monitoring and designing.The aim of this model is to provide a first approximation of interplay of various physicochemical phenomena involved in a RHVFC.It is hoped that this model will accelerate the advancement of the knowledge and quantification of the relevant physical and electrochemical processes, helping to distinguish the key phenomena.The paper presents a detailed explanation of the proposed comprehensive model, stating equations, assumptions, limitations and capabilities of the model to represent experimental data, as such it is expected to provide a baseline for more in-depth representations developed in future.In the following sections the model equations are presented followed by the experimental details and a discussion of the model calibration and validation.

Mathematical Model
The proposed mathematical model was developed considering six domains: (1) catholyte tank, (2) cathode, (3) membrane, (4) CL, (5) Gas Diffusion Layer (GDL), and (6) the anode channel, as shown in Figure 1.The cathode channel was not included as a domain since the flow rate through it was assumed to by-pass the porous electrode and to not participate in any reaction.The hydrogen tank was not considered in the model because the hydrogen flow is passed through the cell and not recirculated to a tank in the experimental set-up that was used (see Experimental section).
The model was developed based on the following assumptions and simplifications: 1.All domains are considered isothermal.2. Electrolyte is considered incompressible, having constant density, viscosity and volume.

F1719
3. Physical properties and mass and charge transfer properties are assumed isotropic and homogeneous in all the domains.4. The electrolytes must maintain electro-neutrality.5. Unit activity coefficients are assumed for all species.However, an activity term is considered for the chemical dissociation of H 2 SO 4 , and a fitted global activity coefficient is considered for the estimation of OCP. 6. Dilute solution approximation is considered.7. Spatial concentration distribution in the catholyte tank, cathode, CL and anode channel are neglected.8. Membrane and GDL are considered to be steady-state.9. Gas evolution reactions in the cathode are neglected.10.Dissociation of bi-sulfate ions is assumed to reach equilibrium instantaneously.11.Full humidification of the membrane is assumed at the cathode -membrane interface.12.The ionomer phase that is present in the CL has the same fixed charge concentration as the ionomer in the membrane.13.Protons are the only ionic species that can cross the membrane.14.Water crossover from cathode to anode is considered, but any change of the catholyte volume is neglected.15.Water only in vapor phase can leave the CL and enters the anode channel or vice-versa, and can leave the system.16.The gas phase at the anode side is treated as an ideal gas.
Cathode and catholyte tank.-Thecatholyte has been considered to contain water, sulphuric acid species (H + , HSO − 4 and SO 2− 4 ) and charged vanadium species (VO 2+ and VO + 2 ).A correct estimation of the initial concentration of ionic species in the catholyte is important to maintain electro-neutrality ( z i c i = 0), and to calculate the equilibrium potential (Equation 40) and overpotential (Equation 46).The H 2 SO 4 dissociation is considered to follow a two-step processes, Equations 6 and 7, where k 1 and k 2 represent the first and second dissociation constants, respectively.
The first step of dissociation is assumed to be fully complete since it has been reported that it is essentially complete for concentrations up to 40 mol kg −1 at temperatures between 273 and 323 K. 57 In chemical equilibrium, the second step of dissociation, that noticeably depends on temperature, is described by the dissociation constant of bi-sulfate ion (k 2 ), which is defined by the activities (a i ) of the sulphuric acid species and can be related to their molar concentrations (c i ) and activity coefficients (γ i ) according to The derivation of Equation 8 is similar to the one presented by Knopf and co-workers, 57 but the reference concentration has been changed from molality to molarity, and includes a molar equilibrium dissociation quotient, Q M (T ), and an activity coefficient product, γ(T ).The experimental molality data of H 2 SO 4 dissociation presented in their work 57 was transformed to molarity and fitted to a simple empirical model.Figure 2 presents the values of Q M at 290 K as a function of H 2 SO 4 molar concentration and the proposed empirical model.
The conservation of mass for each ionic species in the cathode during charge and discharge is presented in Equations 9-13.These equations include the change in the concentration of species entering to the cathode due to the recirculation of the catholyte through the tank; the electron transfer reaction according to Equation 1; the dissociation of bi-sulfate ion according to Equation 7; and the transport of protons through the membrane.The term S eq d that represents the change in concentration of species due to the dissociation of bi-sulfate ions Exp.data model (Equation 7) is obtained by means of Equation 8. Table I shows the input parameters used on the cathode side.
eq d [11] The terms c i and c i,T denote the bulk concentration of species i in the electrolyte at the cathode and catholyte tank, respectively.The symbol ε k , V k and A k represent the porosity, volume and crosssectional are of the domain k.The term t is time, Q ca is the flow rate of catholyte in the cathode, j is the applied current density (current per unit of cross-sectional area of electrode), and F is the Faraday constant.
The mass conservation equations for all ionic species in the catholyte tank is given as follows: Where V T is the catholyte tank volume.
The catholyte flow rate flows through the serpentine flow channel while penetrates into the porous electrode.A steady state condition was assumed to describe a simplified flow distribution inside of the positive half-cell, 59 i.e., p ca = p CH , allowing to calculate the flow rate in the channel and porous electrode.The catholyte flow rate that flows through the flow channels was assumed to by-pass the reaction area (porous electrode).According to Darcy's law, the pressure drop through the porous electrode ( p ca ) can be calculated as, 59 p ca = μL ca Q ca κA ca [15]   where the effective permeability of the cathode (κ) can be determined using Kozeny-Carman equation 30,59 as below, In the above equations, μ is the dynamic viscosity of the catholyte, which was considered constant at its value for 50% State Of Charge (SOC). 30L ca is the cathode flow length, A ca is the flow area of electrode, d f is the mean fiber diameter of the electrode, and K KC is the Kozeny-Carman constant.
The pressure drop through the channel can be calculated considering two sources: friction pressure drop and minor pressure drops, which are caused by sudden change in fluid direction or velocity. 59onsidering these two parts, the pressure drop through the channel ( p CH ) can be obtained as, 59 p where the friction resistance coefficient ( f CH ) can be obtained as shown below, The value 56.91 in Equation 18 is determined considering the specific geometry of the channel. 59,60In the above equations, D h,ch is the hydraulic diameter of the channel (D h,CH = 4A CH /P CH ), K f is the minor resistance coefficient, ρ is the density of the catholyte, A ch is the area of the channel, Q ch is the flow rate in the channel and Re D h is the Reynolds number based on the hydraulic diameter.
Anode catalyst layer.-Theanode side was modeled considering the CL, GDL and anode channel.Figure 3 shows the molar flow rates ( ṅi ) that were considered in the mass conservation equations of the different species.In Figure 3, the direction of the arrows represent positive values for the mass transport of species through the membrane and GDL, reaction of hydrogen in the CL, desorption of liquid water in the CL, and phase change between vapor and liquid water in the CL and anode channel.At the CL, protons ( ṅH + ,m ) and dissolved water can arrive or leave due to transport through the membrane.Significant selfdischarge of redox flow batteries due to crossover of ionic species has been observed, 33,34 which increases with current density.Crossover fluxes of VO 2+ and VO + 2 through the membrane have been neglected in this model (Assumption 13), but will be incorporated in future work.The dissolved water is transported by electro-osmotic drag ( ṅEO ) and diffusion ( ṅdiff ), which was considered to always be in the direction of the CL because of Assumption 11.The dissolved water in the CL can be desorbed ( ṅdes ) into liquid water, which can receive or lose material by evaporation/condensation ( ṅwv,CL ) depending on the vapor pressure and its saturation pressure.Hydrogen/protons are oxidized/reduced by means of the electron-transfer reaction ( ṅH 2 ,react ) and hydrogen can be transported from the CL through the GDL into the channel or vice-versa, along with water vapor ( ṅH 2 ,GDL and ṅv,GDL ).In the channel, pure hydrogen is entering ( ṅH 2 ,IN ), the phase change between vapor and liquid is considered ( ṅwv,CH ), and the outlet only contains hydrogen and water vapor ( ṅH 2 ,OUT and ṅv,OUT ).
The mass conservation equation for each species during charge and discharge in the CL is presented in Equations 19-23.Table II shows the parameter values for the anode side.The concentration of protons (Equation 19) is presented as steady-state because of Assumption 4, which is c H + + z f c f = 0 where z f and c f are the charge number and concentration of the fixed charge in the electrolyte.The balance of hydrogen and water vapor includes the transport through the GDL, which is explained in the following Gas diffusion layer subsection.
In the above equations, V CL is the CL volume, ε m is the volume fraction of electrolyte membrane in the CL, A CL is the cross-sectional  24 area of the CL, ε CL is the porosity of the CL ans s CL is the liquid saturation in the CL.
Water was considered to exist in three different phases depending on the domain, 24,27 i.e., water vapor (v), liquid water (w) and dissolved water (dw).The transport of water through the polymer electrolyte is assumed to happen in the dissolved phase since the pore size of the polymer electrolyte is of the order of only 10 nm and clusters of water tend to be localized and less connected. 24The mass conservation for the water species in the CL are given as follows: Where γ des is the desorption rate, ξ drag is the electro-osmotic drag coefficient, D dw,m is the diffusion coefficient of dissolved water in the electrolyte membrane, l m is the thickness of the membrane, ρ w is the water density, M w is the molar weight of water and R w,CL is the term that represents the phase change between vapor and liquid water in the CL.
A diffusive approach was considered for the transport of dissolved water through the membrane, Equation 21, accounting for electroosmotic drag with ξ drag = 2.5λ/22, 20 diffusion transport and a desorption term.The diffusion coefficient of dissolved water is given by Equation 24.The desorption was assumed to be proportional to the difference between the actual electrolyte membrane water content, λ (Equation 25), and the equilibrium water content, λ eq (Equation 26), which was determined by experimental results of water uptake in CLs. 61The desorption rate was defined by an empirical relationship that depends on the local electrolyte membrane hydration, 24 Equation 27.
In the above equation, ρ dm is the dry membrane density, E W is the equivalent molecular weight of the dry membrane, a i is the activity of species i, f V is the volume fraction of water in the membrane, V w water molar volume, V m is the dry membrane molar volume, l CL is the thickness of the CL, and T is the temperature.V m can be calculated by V m = E W/ρ m . 21he phase change between liquid and vapor for the water is presented in a switching format, Equation 28.Where k c and k e are the condensation and evaporation rate constants, y v is the molar fraction of vapor in the gas phase, p v is the vapor pressure and p sat v is the saturation pressure.The mass-transfer rate is related to the amount of reactant in the porous media and the driving force, 25 which is the difference between the vapor pressure and its saturation pressure.The restriction that when no liquid water is present and the vapor pressure is lower than the saturation pressure, Equation 28is not valid and no evaporation can occur, must be included.
) [28]   Gas diffusion layer.-TheDusty Gas Model (DGM), which considers both diffusion and pressure-driven convection, 64 was used to model transport through the GDL and only the steady-state transport of hydrogen and water vapor was modeled.Equation 29 was used to describe the transport of a binary mixture 65 (H 2 and water vapor (v)), accounting for the difference of molar fraction and pressure.The molar flow rate of species can be obtained by means of the following relation, ṅi,GDL = A GDL N i,GDL , where A GDL is the cross-sectional area of the GDL and N i is the molar flux of species i.The effective binary diffusion coefficient (D eff i,j ) was estimated from Chapman-Enskog kinetic theory assuming ideal gas behavior holds 66 (Assumption 16) and is given in equation 30.The effective Knudsen diffusion coefficient (D eff Kn,i ) and permeability of the porous material (K ) were calculated using Equation 31. 65 Both diffusion coefficients were corrected considering porosity and tortuosity (τ) of the porous media using a Bruggeman correlation, 19,67  where, Anode channel.-At the inlet, dry hydrogen is entering, and at the outlet, hydrogen and water vapor leave the anode channel.Liquid water was not considered to leave the channel (Assumption 15).The mass conservation equations of hydrogen, vapor water and liquid water are given in Equations 35, 36 and 37, respectively.Inside the anode channel the phase change between vapor and liquid water is given by Equation 28.The mass flux (G) that enters and leaves the anode channel was estimated through an isothermal energy balance, 69 Equation 38, considering the fanning friction factor ( f F ) for laminar flow, f F = 16/Re, where Re is the Reynolds number.Figure 3 shows the two positions considered in each energy balance in the inlet ( p 1 and p 2 ) and outlet ( p 1 and p 2 ).After solving Equation 38 assuming laminar flow at the inlet and outlet, the Reynolds number was recalculated and a value around 4.5 was obtained for inlet and outlet, which fulfills the previous assumption (i.e., Re < 2000).
In the above equations, V CH is the anode channel volume, s CH is the liquid saturation in the anode channel, R w,CH is the term that represents the phase change between vapor and liquid water in the anode channel, M is the molecular weight of the gas, d is the pipe diameter, L is the pipe length, p 1 and p 2 are the pressure and the initial and final point, ρ is the gas density, and v m is the mean velocity of the fluid.
Cell potential.-Theoperating cell potential, E cell was estimated considering the reversible OCP (E OCP ), ohmic overpotential (η ohm ) and electrode overpotentials (η ca ) and (η an ), 70 Equation 39.Where the overpotentials are added to E OCP for charge operation and are subtracted for discharge operation.
A complete Nernst equation to estimate the OCP was derived form thermodynamics as explained in Appendix A. Considering Equation A12 and assuming the activity of water (a ca H 2 O ) in the catholyte equal to 1 and Henry's law for hydrogen, which considers the activity of dissolved hydrogen equivalent to the change in partial pressure of hydrogen. 35The OCP can be expressed as, To further simplify Equation 40 and considering the difficulty in obtaining data of activity coefficients, a global activity term was considered (γ OCP = γ ca ). Equation 40 considers the potential difference of the catholyte and anolyte at equilibrium, which was obtained by the condition of equality of electrochemical potential of protons in the catholyte and anolyte, [42][43][44] see Appendix A. An equivalent expression for OCP can be derived, if the Nernst equation is considered along with the potential jump (Donnan potential) that appears at the two membrane -electrolyte interfaces due to the difference in concentration of protons. 42This potential across the entire membrane, namely dialysis potential, derived from thermodynamics differs from the relation presented by Knehr and Kumbur 38 as has been stated by Pavelka et al. 42 The ohmic overpotential was calculated by adding the contribution of the different domain: current collectors (cc), catholyte (ca-el), cathode (ca), membrane (m), CL, polymer electrolyte in the CL (CL-m) and GDL; and were described by means of Ohm's law, 46,72 as shown in Equations 41 and 42.In addition, an extra ohmic resistance (R C ) was included to account for any additional resistance, such as the contact resistance between current collectors and porous electrodes 10,73 and changes in the membrane conductivity due to uptake of ionic Where l k and σ k are the thickness and conductivity of the domain k, and the superscript 'eff' denotes an effective property.Table IV shows the property values for the membrane and current collectors.
The Nafion membrane conductivity was described as a function of the water content (λ) and temperature through an empirical relation presented in Equation 43. 20,72m = (0.5139λ − 0.326)exp 1268 The effect of the membrane uptake of sulphuric acid and vanadyl on membrane conductivity has been studied by Tang et al., 74 who observed an enhancement or reduction of membrane conductivity depending on acid and water content in the membrane.An increase in the concentration of protons due to acid presence in the membrane can enhance conductivity, while mobility of protons significantly decreases with water loss. 74When a nafion membrane was equilibrated with a electrolyte solution of practical composition for VRFB operation, uptake of vanadyl and acid contributes to reduce membrane conductivity to some extent, and a complicated equilibrium and partitioning competition between protons and vanadyl has been suggested. 74Future studies should focus on the membrane-electrolyte equilibrium in the RHVFC to assess membrane conductivity, concentration and mobility of protons by means of a detailed transport model of water and ionic species through the membrane.
The electrolyte conductivity was estimated considering the concentration and charge of all ionic species in the electrolyte, 45 Equation 44.A Bruggemann correction was applied to account for the porous electrode, 45 allowing for the calculation of the effective conductivity of electrolytes and electrodes as shown in Equation 45.
In the above equations, z i is the charge number of species i and D i,el is the diffusion coefficient of species i in the electrolyte.
The electrode overpotential with contribution of cathode and anode was considered to estimate the cell potential.The cathode overpotential (η ca ) was calculated by a Butler-Volmer (BV) kinetic approach 45,46 (Equation 46), where all the ionic species involved in the redox reaction were considered. 75,76This expression does not consider activity coefficients for the ionic species, because of the difficulty in estimating individual activity coefficients or decoupling their values from the fitted global activity terms (γ OCP ).It is expected that the fitting parameters used in the Butler-Volmer equation absorb to some extend the effect of neglecting the activity coefficients.
The applied current density can be obtained according to j = S ca V ca j BV /A ca , and the exchange current density was estimated by Equation 47.The dependency of the cathode rate constant with temperature can be expressed by an Arrhenius approach, Equation 48. [47] In the above equations, the term j BV and j BV 0,ca are the current density and exchange current density (current per unit of surface area of pore walls), α is the transfer coefficient for the vanadium redox reaction (Equation 1), S ca is the specific surface area of the cathode (surface area of pore walls per unit volume of total electrode) and c s i is the surface concentration of species i at the liquid-solid interface of the cathode.k ca and k ca,ref are the rate constant and reference rate constant, T ref is the reference temperature and E 0 ca,T ref is the standard potential at the reference temperature.
The Butler-Volmer equation presented in Equations 46 and 47 considers the mass transfer effects, which become more important at higher applied current densities or even for a small current density at high SOC conditions during charge-discharge operation.In order to evaluate the surface concentrations, the rate of electrode reaction (i.e., flux of species consumption/production according to Butler-Volmer equation) was matched to the rate at which electro-active species are brought to the surface by mass transfer, 76 Equations 49, 50 and 51.The rate of mass transfer was considered proportional to a linear concentration gradient at the surface of the electrode within a Nernst diffusion layer (δ 0 ). 30,41,45,76

F1724
Journal of The Electrochemical Society, 164 (14) F1717-F1732 (2017) The overpotential for the anode (η an ) was approximated by a Tafel-Volmer approach, which was utilized to describe the hydrogen oxidation/evolution reaction, since it has been reported that these reactions are not successfully represented by a Butler-Volmer relation. 35Kinetic models considering two pathways with three elementary steps (i.e., Tafel, Heyrovsky and Volmer) are well established to develop the overall rate expression.Tafel-Volmer (TV) elementary steps include adsorption/desorption of reactant and product and the electron transfer reaction, 35 which are shown in Equations 52 and 53.
Where S is the surface adsorption site, H 2,surf is the hydrogen molecule adjacent to the surface, k ad and k des are the adsorption and desorption rate constants of Tafel reaction, and k V and k −V are the forward and backward rate constants of the Volmer reaction.
A TV kinetic model that has been introduced by Kucernak et al. 35 was utilized in the proposed model, Equation 54.This equation describes the current density -overpotential relation as independent of the pH, and only dependent on the hydrogen partial pressure and kinetic rate constants under near mass-transport free conditions.The coverage of hydrogen on the electrode surface (θ TV H ad ) is also a function of the hydrogen pressure and the overpotential (Equation 55).Crossover of vanadium species is expected to reduce the electroactive area in the anode by means of physisorption on the Pt catalyst competing with hydrogen adsorption, which will affect the current density of the cell.This loss in electro-active area must be included in the model by modifying the coverage of hydrogen on the electrode surface (θ TV H ad ) and/or the roughness factor of the anode (R an ); however such an effect is not considered in the present model.The effect of temperature was considered in Equation 57, which assumes that the reaction is controlled by a single activation energy (E a ). 35he parameters of the TV model reported at standard ambient temperature and pressure (SATP, ) conditions, namely a temperature of 298.15K and an absolute pressure of 1 bar, are given in Table III.The current density at standard conditions ( j TV, ) is obtained using Equations 54 to 55 and the applied current density can be obtained according to j = R an j TV , where R an is the roughness factor of the CL and j TV is the current density at operating conditions (current per unit of electrochemically active area).
where, [57]   In the above equations, β is the transfer coefficient of the hydrogen redox reaction (Equation 2) and a H 2 is the activity of dissolved hydrogen and equivalent to the change in partial pressure of hydrogen 35 according to a H 2 = a H 2 ( p H 2 / p H 2 ), where a H 2 and p H 2 are the standard state activity of hydrogen in solution and in the gas phase (a H 2 = 1, p H 2 = 1 bar).

Experimental
Experimental data was collected to validate the proposed model, galvanostatic tests were performed using an in-house manufactured RHVFC with a cross-sectional area of 25 cm 2 .The cell consisted of polypropylene insulating layers, aluminum end plates, graphite flow channel plates and the MEA.The MEA contained a carbon electrode (SGL 10AA, 400 μm) and a platinised carbon paper electrode (Alfa Aesar Hydrogen Electrode, 200 μm, 0.5 mg cm −2 Pt loading), which were separated by a Nafion 117 membrane.Multi-channel serpentine flow channel plates (SGL, BMA5 graphite plates) were used as current collectors and to distribute the vanadium electrolyte as well as the hydrogen gas.To improve the current collection, copper current collectors were used between the polypropylene plates and flow channel plates in a previous set of experimental tests. 77The catholyte solution was prepared by fully dissolving 14.1 g of vanadium sulfate hydrate (Sigma-Aldrich) in 60 mL of 5 M H 2 SO 4 solution (Fluka Analytical).A peristaltic pump was used to circulate the vanadium catholyte between the cell and catholyte tank at a constant flow rate of 1.67 × 10 −6 m 3 s −1 .Hydrogen was passed through the anode side at a constant flow rate of 1.67 × 10 −6 m 3 s −1 .A collector reservoir was connected to the hydrogen outlet to collect any catholyte crossover.
The galvanostatic charge and discharge tests were performed using a Bio-Logic potentiostat (VSP-300) running EC-Lab software.Experiments were performed at constant current density in the range of 50-600 A m −2 and flow rate of catholyte and hydrogen in the range of 0.42-2.5 × 10 −6 m 3 s −1 , allowing the cell to reach an upper cutoff potential of 1.3 V and a lower cutoff potential of 0.4 V. OCP measurements were performed after each galvanostatic charge-discharge step.Finally, the OCP behavior as a function of SOC was measured by charging or discharging on a series of steps at constant current density and measuring OCP after each step.The SOC of the cell was calculated by comparing the experimental capacity with the maximum theoretical capacity (100% SOC) and considering that fresh solution was used ( 0% SOC).A summary of the operating conditions used in the OCP test and the charge-discharge (ch-dch) tests used in this paper is presented in Table V.

Results and Discussion
The model equations described in Mathematical model section were developed and solved in MATLAB by means of an Ordinary Differential Equation (ODE) solver, with absolute and relative tolerance set at 1 × 10 −6 .Additionally, the events option of the ODE solver was used to incorporate the restriction associated to Equation 28. 78The simulation of a potential curve in Figure 5   catholyte, 38 which were estimated as follows: Where c V is the total concentration of vanadium and c 0 i is the initial concentration of species i.
The proton, bisulphate and sulfate concentrations were assumed to follow a complete first step of dissociation of sulphuric acid (Equation 6), and a second step (Equation 7) according to a bi-sulfate dissociation constant (k 2 ) and an activity coefficient product (γ) as presented in Equation 8.For this purpose, the molar equilibrium dissociation quotient, Q M = k 2 /γ presented in Figure 2 was used.
Figure 4a shows the E OCP estimated by Equation 40 compared to one data set of experimental OCP.This Figure includes three different estimations: (1) a the Nernst Equation (NE), Equation 63, with unity activity coefficients, (2) a Complete Nernst Equation (CNE), Equation 40, with unity activity coefficients, and (3) a CNE considering non-ideal activity coefficients by means of a global activity coefficient (γ OCP ).The fitted values of γ OCP against experimental data are shown is Figure 4b.
In the case of VRFBs, the Nernst equation with unity activity coefficients underestimates the OCP data and the Donnan potential across the membrane has been added in order to fit OCP data. 38However, this potential difference should reduce the OCP when the catholyte concentration of protons is higher than the anolyte one, because a steady state condition must be maintained, i.e., no net flux of protons across the membrane. 42This effect is contrary to the one presented by Knehr and Kumbur, 38 which increased the OCP in the case of all-vanadium cells because the equation considers the catholyte concentration of protons to the power of 3. On the other hand, Pavelka et al. 42 have presented a thermodynamic derivation of OCP for VRFBs that does not contain this cubic term, but instead the catholyte concentration of protons to the power of 1.They stated that for cation exchange membranes the non-ideal activity coefficients should be found or fitted against experimental data to improve the agreement. 42he OCP equation used in this work was derived as shown in Appendix A and is analogous to the derivation presented by Pavelka et al. 42 Similar effect to the described above is shown in Figure 4a, where the CNE with unity activity coefficients gives lower values of OCP in comparison to the NE curve.The fitted values of γ OCP are shown in Figure 4b along with an empirical model of its dependency  with SOC.The fitted values of the global activity coefficient vary between 10 to which be related to the value of the activity coefficient of protons that have been reported to be of approximately 5 for a concentration of H 2 SO 4 of 1 mol kg −1 and even higher for higher concentrations of H 2 SO 4 . 57This analysis was repeated for each OCP experimental data set, which showed similar values for the fitted global activity coefficient.OCP experimental measurements were performed four times considering two different fresh solutions.A correct estimation of the OCP is key for the accuracy of model predictions and the estimation of overpotential contributions.
Model validation -potential.-Thesecond step in the model validation is to compare the model prediction against the experimental potential data of charge and discharge.The fitting procedure was implemented in MATLAB using a non-linear least-squares solver (lsqcurvefit function), which allowed upper and lower boundaries for the fitting parameters to be set.In principle any model parameter could be taken as a free fitting parameter, but its value should be maintained within the physically meaningful limits. 79Common fitting parameters used in VRFB models are reaction rate constants, ionic or electronic conductivities, electrode specific surface areas, charge transfer coefficients, diffusion coefficients and permeabilities. 10,30,39,45,46It is important to select realistic initial values of the fitting parameters to decrease the solution uncertainty since the solver algorithm does not necessarily find a global optimum or unique solution. 79To this end, the fitting parameters were chosen to be the standard reaction rate constant for the cathode (k ca,ref ), the Nernst diffusion layer (δ 0 ) and an extra ohmic resistance (R C ). Their initial values were selected as k ca,ref =8.1 × 10 −7 ms −1 , 80 δ 0 =50.3 × 10 −6 m which is the mean pore radius (r p ), 30 and R C =5 cm 27 for the cell without Cu current collectors and R C =0.5 cm 2 for the cell with Cu current collectors.Their lower and upper boundaries were selected as 10 −10 and 10 −6 m s −1 for k ca,ref , r p /10 and 2r p for δ 0 and, 3 and 6 cm 2 for R C for the cell without Cu current collectors, while 0.01 and 1 cm 2 for R C for the cell with Cu current collectors.The decision of including an additional resistance was based in previous measurements 7 which showed a series resistance (R S ) of around 5 cm 2 at a catholyte flow rate of 1.67×10 −6 m 3 s −1 using the same experimental set-up (cell, assembly, supplies, etc.) that the one used in this work without Cu current collectors, while a series resistance of around 0.8 cm 2 with Cu current collectors.
The fitting of the model was carried out against one single cycle experimental data at a current density of 400 A m −2 and then simulations were carried out using the obtained fitting parameters to validate the model under different current densities, as reported in Figures 5 to 7. The only parameter refitted was R C , depending if the tests were done with or without Cu current collectors.The simulations were performed following the actual operation of the cell, setting the initial SOC at a very small value (0.1%) in the cases where fresh catholyte solution was used or to an estimate value of initial SOC usually between 0 and 5%, calculated based on the experimental data.
The initial concentration of species and the operating parameters that were considered in the simulation are given in Table VI.
Figure 5 presents the model calibration against experimental charge and discharge potential curves at a current density of 400 A m −2 , Table V data set N • 5, for the cell with Cu current collectors.During the model fitting, it was noticed that the fitting error was mainly due to the differences at the final part of the discharge curve (SOC<20%), and therefore, larger weighting factors were attributed to the charge curve and the beginning of the discharge curve during the model fitting process to allow for a more sensible identification of fitting parameters.The fitted parameters obtained were k ca,ref =1.2 × 10 −7 ms −1 , δ 0 =84.8 × 10 −6 m and R C =0.3 cm 2 .Although, the best fit could not remove all the discrepancies between the experimental data and the model simulation, a reasonably good agreement was observed with a total RMSE of 3.6% and 5.4% during charge and discharge, respectively.It can be seen that the error for the charge simulation is mainly produced by differences at extreme conditions and the high error for the discharge simulation is mainly due to the differences at low SOC, which are less critical for practical applications.It is worth noting that the rate of depletion of active species in the cathode is sensitive to the flow rate penetration into the cathode, which was calculated by a simplified model (Mathematical model section) along with the assumption that only the fraction of catholyte flow rate through the cathode participate in the reaction.This may well be responsible for most of the differences during discharge.Additionally, the assumption of unity activity coefficients (Assumption 5) and neglected effects such as the transport of vanadium and H 2 SO 4 species through the membrane into the CL and changes of water content in the catholyte could explain to some degree the difference between the model and experimental data.Self discharge and gas evolving reactions are not considered in the model, which are expected to decrease the coulombic efficiency in the experimental case.
After the model calibration, the parameters k ca,ref , δ 0 and R C were used to simulate the cell performance at 80 and 600 A m −2 , Table V data set N • 4 and N • 8, respectively.The model demonstrated good accuracy in representing the performance of the cell at both current densities.Figure 6 shows the validation of the model by comparing the simulated cell performance and the experimental data at 8 and 600 A m −2 , which produced a RMSE of 0.3% and 3.8% during charge and discharge at 8 A m −2 and a RMSE of 1.4% and 3.8% during charge and discharge at 600 A m −2 .An increase in the discrepancy between the model and the experimental data occurs at the beginning of discharge when a higher current density was used, which may respond to an overestimation of the average concentration of species in the cathode.The parameter R C was refitted to be able to simulate the experimental data at 50 and 100 A m −2 of the cell without Cu current collector, Table VI data sets N • 2 and N • 3. The refitted R C took the value R C =3.6 cm 2 .The model demonstrated good agreement with the experimental data at both current densities.Figure 7 shows the validation of the model by comparing the simulated cell performance and the experimental data at 50 and 100 A m −2 , which produced a RMSE of 1.6% and 2.0% during charge and discharge at 50 A m −2 and a RMSE of 0.7% and 1.8% during charge and discharge at 100 A m −2 .It can be noticed the effect of the extra ohmic resistance, which produce an increase of approximately 60 mV of the potential curve at 100 A m −2 compared with the curve at 80 A m −2 over the linear range of performance.On the other hand, validation of the model at different flow rates of catholyte and hydrogen (data sets N • 6 and N • 7) can be seen in Figure 8.The model was able to simulate the performance of the cell when different flow rates were used.Evolution of species in electrodes.-Theresults of the evolution of species in the cathode and anode for the experimental data at a current density of 50 A m −2 (data set N • 2) is presented in Figures 9 and  10.The evolution of the species concentration in the cathode during charge and discharge mode is presented in Figure 9.The concentration of vanadium species VO 2+ and VO + 2 based on Equations 9 and 10 and the concentration of sulphuric species based on Equations 11, 12 and 13 are presented in Figures 9a and 9b, respectively.The species concentrations presented a linear rate of change, which is reached at around 70 s for the operating conditions.The rate of change of the species concentrations depends on the stoichiometric coefficients and number of electrons consider in the redox reactions, along with the flow rate penetration to the domain as it can be noticed in Equation 9to 13.In the particular case of protons, 2 moles of H + are generated in the catholyte with respect to 1 mole of VO + 2 during charge, and it would be expected that its concentration increases with this same proportion.However, the concentration of protons in the catholyte is affected because 1 mole of protons must be transported across the membrane to participate in the anode reaction and the chemical dis-sociation of HSO − 4 also takes place in the electrolyte.It can be seen in Figure 9b that the concentration of SO 2− 4 and HSO − 4 change during charge and discharge operation, which is explained by the condition of equilibrium of the reaction of HSO − 4 dissociation, Equation 7. On the other hand, the effect of the flow rate penetration can be illustrated if a lower flow rate penetration is allowed, a faster rate of change of ionic concentrations would be expected.The error produced by the simplified model used to obtain the flow rate penetration into the cathode is acceptable for the discrepancies observed between experimental data sets and the simulations (Figures 5 to 8).The error is expected to increase if higher current densities are used, in which case the model should be extended to include spatial distribution effects, i.e., One or Two Dimensional (1D or 2D) model, allowing the actual species flux into the cathode and across the membrane to be calculated.
The evolution of the species concentration in the anode CL during charge and discharge mode at a current density of 50 A m −2 is presented in Figure 10 and 10b, respectively.All species concentrations in the anode CL presented a steady-state behavior, which is reached in around 10 s for the operating conditions.This initial transient behavior is in response to the water transport through the membrane and the phase change processes of liquid water desorption and evaporation/condensation.It was observed that during discharge the liquid saturation slightly decreased with respect to its value during charge, possibly due to the change of direction of the electro-osmotic drag flux toward the cathode while the diffusion flux maintained its direction toward the anode.It is worth noting that the simulation results somewhat overestimated the amount of water lost from the catholyte tank with respect to the experimental observations.During the experimental operation of a single cycle of charge and discharge at a low current density, the total crossover solution collected from the outlet of the anode channel was marginal compared with the catholyte volume.The dark blue color of this crossover solution suggested that vanadium species were also transported across the anode side and left the system.Tang et al. 74 have observed sulphuric acid and vanadium uptake in membranes after equilibration with solutions of sulphuric acid and vanadyl sulfate.They reported a reduction to some extend of the membrane conductivity, proton mobility and water content in electrolytes with a typical composition of VRFB feed, suggesting a complicated equilibrium and partitioning competition between protons and vanadyl. 74A detailed model of the simultaneous transport of species across the membrane was not explored in this study.Moreover, it is important to note that if different operating conditions were used as higher current densities, it is possible to produce flooding of the anode side.Further study is currently ongoing to describe the water and ionic species transport across the membrane and to estimate their effect on the anode performance.

Conclusions
This study introduced a time-dependent model for a RHVFC considering mass conservation and electrochemical processes.A wellestablished modeling approach has been used to describe a RHVFC by means of coupling physicochemical phenomena used to simulate the performance of different systems, such as VRFBs and PEMFCs.The model was validated considering experimental data of OCP and cell potential.The OCP data was fitted with a complete Nernst equation that was derived from thermodynamics, considering a global activity coefficient.The simulated cell potential, considering the overpotential of cathode and anode and the ohmic losses, was compared to charge-discharge cycle data sets showing a good degree of accuracy in predicting the cell performance.The discrepancy between experimental data and the model simulations at the end of charge and discharge is most likely explained by the the use of a simplified model to estimate the catholyte flow rate penetration into the cathode domain, which may predict a slower rate of change of active species concentration.The discrepancies at extreme conditions (10%>SOC or SOC>90%) could be explained by the use of unity activity coefficients in the kinetic relations, i.e., Butler-Volmer relation for the positive side and Tafel-Volmer relation for the negative side.Moreover, self discharge, gas evolving reactions and crossover of ionic species through the membrane have been neglected in the proposed model, which may well affect the performance of the cell.The model is capable of representing the voltage dynamics observed in a RHVFC at moderate current densities, considering the range of operating conditions used to validate it.However, the use of the model to predict performance at conditions beyond the validation range implemented in this paper might not produce meaningful results and significant further experimental data would be required.
The model presented in this study can be used as a first approximation, allowing simulation of the system and providing a foundation for further development of physical-based models for regenerative fuel cells.In addition, the electrochemical approaches used in the model may serve as reference for studying similar systems, such as allvanadium and H 2 /Br 2 , which may benefit from the complete Butler-Volmer equation for the cathode considering VO 2+ , VO + 2 and H + , the complete Nernst equation and the kinetic approach for the hydrogen oxidation/evolution reaction used in this study.Further research will include RHVFC testing under a wider range of operating conditions and the corresponding model validation.The contribution of effects such as spatial distribution in the electrodes, ionic species crossover, and water transport through the anode, will be significant under high current density operation.Modifications to the proposed model by incorporating a more detailed water management model, transport of cathode ionic species through the membrane and their interaction with the CL as well as study on species spatial distribution (1D to 2D approaches) need to be further investigated.The low-complexity modeling approach used in this paper has enhanced the understanding of the system performance by coupling physical and electrochemical processes occurring in the RHVFC, enabling identification of the key phenomena and highlighting areas requiring future in-depth study.and EP/K002252/1 for funding this work.We would like to express our appreciation to Dr. Antonio Bertei for all the useful discussions and constant support during the course of this research, more specifically his guidance in the derivation of the electrochemical relations and his comments that helped to improve the manuscript significantly.We are also immensely grateful to Dr. Samuel Cooper for his comments on an early version of the manuscript.

Appendix A: OCP Derivation
The measured OCP is the difference in electrochemical potential of electrons in the two terminals of a battery. 42,81The electrochemical potential of a species i ( μi ) is given by Equation A1, which considers the chemical potential of species i (μ i ) and the effect of potential (φ) on a charged species.
Where the chemical potential of species i can be expressed in terms of the standard chemical potential (μ • i ) and activity of species i as, At equilibrium, each half-cell reaction presented in Equations 1 and 2 can be written in terms of electrochemical potentials as follows: Considering, the equilibrium condition between the phases in contact and using Equation A1 to express the electrochemical potentials in the above equations.Where M and M represent the phase of the wires used to connect the voltmeter to the terminals of the battery.Subtracting Equation A7 to Equation A6.
The difference of potential of the electrolytes (Equation A9) can be obtained recognizing the steady state condition, 42 which leads to the equality of electrochemical potential of protons between both electrolytes.
This same relation can be obtained by considering the Donnan potential in the two electrolyte -membrane interfaces. 43,44At equilibrium, the electrochemical potential of protons in the electrolyte and membrane must be the same.This was expressed for each interface as, To write the cell OCP, Equation A9 was substituted in Equation A8, and the chemical potentials were expressed with respect to Equation A2.where the standard cell potential was defined as, 42 ORCID H. Hewa Dewage https://orcid.org/0000-0002-0089-103X took about 4 min on an Intel Xeon E5-1620v3, 64-bit workstation with 32 GB RAM.In the following subsections: Model validation -open circuit potential and Model validation -cell potential, the complete Nernst equation is fitted to experimental data of OCP, and the proposed model is validated against experimental data of cell potential during charge and discharge.Model validation -open circuit potential.-Thefirst step in the model validation is to fit the proposed complete Nernst equation, Equation 40, to the measured OCP.The initial species concentrations in the catholyte were calculated considering the vanadium sulfate hydrate, VOSO 4 • xH 2 O, contained 3.5 molecules of water.During operation, the electron-transfer reaction and acid dissociation cause a change in the concentration of all species in the

Figure 4 .
Figure 4. OCP at operating conditions presented in Table VI.(a) shows a comparison of Nernst equation (Equation63) with unity activity coefficients, and a complete Nernst equation (Equation40) with unity activity coefficients and fitted values of a global activity coefficient (γ OCP ), and (b) shows the fitted values of the global activity coefficient.

Figure 5 .
Figure 5. Model calibration against experimental data at a current density of 400 A m −2 and flow rate of catholyte and hydrogen of 1.67 × 10 −6 m 3 s −1 , data set N • 5.

Figure 6 .
Figure 6.Model validation by using the same fitted parameters against experimental data at a flow rate of catholyte and hydrogen of 1.67 × 10 −6 m 3 s −1 and a current density of (a) 80 A m −2 data set N • 4, and (b) 600 A m −2 data set N • 8.

Figure 7 .Figure 8 .
Figure 7. Model validation, using refitted R C , against experimental data at a flow rate of catholyte and hydrogen of 1.67 × 10 −6 m 3 s −1 and a current density of (a) 50 A m −2 data set N • 2, and (b) 100 A m −2 data set N • 3.

Figure 9 .Figure 10 .
Figure 9. Evolution of species concentration against time in the cathode during charge and discharge at a current density of 50 A m −2 and a flow rate of catholyte and hydrogen of 1.67 × 10 −6 m 3 s −1 , data set N • 2: (a) shows the vanadium species VO 2+ and VO + 2 and (b) shows the sulphuric acid species H + , SO 2− 4 and HSO − 4 .

Table III . Electrochemical and transport parameters.
−1/2.The effect of liquid saturation over the transport parameters was neglected because of Assumption 15, which results in no liquid water in the GDL.In the above equations, R is the universal gas constant, K is the permeability, μ is the gas dynamic viscosity, M i is the molecular weight of species i, σ ij is the characteristic binary Lennard-Jones length, D is the diffusion collision integral, and d p is the mean pore size.The parameters σ ij and D are given by Equations 32 and 33, where σ i and i are the characteristic Lennard-Jones energy and length, ij is the characteristic binary Lennard-Jones energy, and k B is the Boltzmann constant.TableIIIpresents the electrochemical and transport parameters that were considered in the model.

Table V . Experimental data sets.
Experimental test: OCP refers to open circuit potential test; and ch-dch refers to a test of a single charge-discharge cycle. a

Table VI . Operating parameters and initial conditions.
, mol cm −2 s −1 k des desorption rate constant, mol cm −2 s −1 k V forward rate constant of Volmer reaction, mol cm −2 s −1 k −V backward rate constant of Volmer reaction, mol cm −2 s −1 R s series resistance, cm 2 Runiversal gas constant, 8.314 J mol −1 K −1 R an roughness factor of CL, m 2 m −2 w phase change rate, mol s −1