A Model for Analysis of the Porous Nickel Electrode Polarization in the Molten Carbonate Electrolysis Cell

It is important to know the electrode kinetics of hydrogen production as well as to understand the effect of the mass transport in the gas phase for the analysis of the molten carbonate electrolysi ...

Because of the limited resources of fossil fuels, and the environmental issues involved with combustion of these, there is currently a growing focus on the production of hydrogen and syngas as alternative energy carriers. 1 A promising option for producing hydrogen and/or syngas is by water electrolysis and/or co-electrolysis of water and carbon dioxide using the molten carbonate electrolysis cell (MCEC). Additionally, with the increased share of power generated from renewable energy sources such as wind and solar, which is fluctuating and intermittent by nature, there is also a need for new energy storage technologies. 2,3 To use the surplus from renewable electricity for the MCEC to produce hydrogen and syngas, which can be stored, is therefore a very attractive possibility.
The molten carbonate electrolysis cell operates at high temperatures, 600-700 • C, which gives the potential to attain high overall efficiency and to require low applied voltage. Our previous study 4 proved that it is feasible to operate the molten carbonate electrolysis cell using the corresponding fuel cell (MCFC) catalysts, i.e., nickelbased porous electrodes. Moreover, lower polarization losses were found for the electrolysis cell than for the fuel cell, 4 and the electrolysis cell also showed good durability in the long-term test using a conventional fuel cell setup. 5 The molten carbonate fuel cell is already commercially available. The possibility of operating both the electrolysis cell for fuel gas production and the fuel cell for electricity generation in a single device may be beneficial from an economic perspective.
In the MCEC, the electrolysis of water takes place on the Ni electrode according to [1] Due to the presence of carbon dioxide in the inlet gas, direct CO 2 electrolysis generating CO may take place at the Ni electrode through the following reaction 2CO 2 + 2e − CO + CO 2− 3 [2] However, the electrolysis of carbon dioxide for CO production is much slower than the electrolysis of water for H 2 production on the nickel-based electrodes. 6,7 The possibility of direct carbon dioxide electrolysis producing CO is small, at least at low current densities. * Electrochemical Society Member. z E-mail: lanhu@kth.se Despite that, carbon monoxide could be generated through the following water gas shift reaction (WGSR).
It has been reported that the shift reaction could rapidly establish equilibrium in the presence of suitable catalysts such as nickel at high temperatures, above 500 • C. 8,9 The presence of the WGSR would affect the gas composition in the electrode and therefore influence the electrode performance.
In order to predict the performance of the Ni electrode under various operating conditions in the molten carbonate electrolysis cell, it is important to know the electrode kinetics of hydrogen production as well as to understand the effect of the mass transport in the gas and electrolyte phases. In the previous study, 10 the experimental results on partial pressure dependencies of exchange current densities have shown higher values than the theoretical dependence obtained from the proposed reaction mechanism of hydrogen production. However, the proposed reaction mechanisms concern planar electrodes, which are different from porous electrodes when regarding electrode surface and mass-transfer limitations. No mechanism has been found to fully explain hydrogen production on the porous electrode in the MCEC.
When determining the kinetics of the hydrogen production reaction in the MCEC, the porous Ni electrode was assumed to be mainly under kinetic control and the influence of mass-transfer limitations in the gas phase was not taken into account. 10 Additionally, the WGSR was found to have an effect on the water partial pressure dependency. 10 The inlet gases used for the Ni electrode in electrolysis cell operation contain H 2 O, CO 2 , H 2 and also N 2 , therefore it is necessary to consider the WGSR.
In this paper, we developed a mathematical model for the MCEC porous Ni electrode based on the Maxwell-Stefan formulation for mass transfer. The model is combined with equations describing the current distribution in the electrolyte phase in the porous Ni electrode. Some previously obtained experimental data 10 for the Ni electrode are used as input for the model. The purpose of the study is to attain knowledge on limiting phenomena, achieve a model able to predict cell performance for a wide parameter space, and to investigate the importance and impact of the WGSR in the cell.

Model
The model is one-dimensional for an isothermal and isobaric system. The model includes the gas-phase mass transport in the current collector and in the porous electrode. When including WGSR, the shift reaction is assumed to mainly take place in the electrode and to be negligible in the current collector, due to much larger specific surface area inside the electrode than in the current collector. In this model the mass transport in the electrolyte phase is not taken into account. The mass transport of species i in the gas phase can be described by the Maxwell-Stefan equation. 11 where R is the gas constant, T is the cell temperature, P is the pressure, x i is the molar fraction of species i and the sum of molar fractions is one, x i = 1, N i is the molar flux of species i with respect to the stationary axis, and D eff i j is the effective binary diffusion coefficient between gas species i and j. Effective diffusion coefficients, following the expression D eff i j = ε τ D i j , were used to account for the effect of porosity and tortuosity in the current collector and porous electrode. For the current collector the tortuosity was set to 1, and for the porous electrode the expression τ = ε −0.5 was used. The binary diffusion coefficient values D i j were calculated using the Chapman-Enskog theory. 11 For the diffusion coefficients involving polar molecules, for which the Chapman-Enskog may be less accurate, the values where compared to corresponding experimental values. 12 Generally good agreement was found, except for the case of H 2 O-CO 2 , where the reported experimental value was about 75% larger than the calculated one. For H 2 O-CO 2 the experimental value was used, whereas the Chapman-Enskog values were used in all other cases. The values used in the simulations are shown in Table I.
The molar flux, including the diffusion and convection, is given by where M i is the molar mass of species i and u is the velocity. A mass balance for gas species i in the electrode is given by where R i is the molar rate of production of i per unit volume. When both the electrochemical reaction of hydrogen production and the water gas shift reaction take place in the electrode, the reaction rate is given by [7] In this study the electrode is regarded as a solid porous matrix where S is the internal specific surface area. 10,13 n is the number of electrons in the electrochemical reaction, F is Faraday's constant, j loc is the local current density from the hydrogen production reaction, r shift is the local reaction rate for the water gas shift reaction, and υ i is the stoichiometric coefficient of species i in the hydrogen production and the shift reaction.
In this study the inlet gases containing H 2 O, CO 2 , H 2 and N 2 are used for calculations (gas compositions listed in Table II). At steady state the time-dependent concentration change of any gas species is zero, i.e. ∂c i /∂t = 0, resulting in the following main equations for the electrode. [10] ∇ · N CO = r shift [11] Furthermore, at steady state there is no flux of nitrogen, since it is inert and is neither produced nor consumed in the electrode.
When considering the water gas shift reaction the reaction rate variable r shift is solved for numerically in the electrode domain as a dependent variable in the model in order to fulfill the equilibrium Here the equilibrium constant K p has been calculated in the previous study. 10 When not including the WGSR, r shift is simply set to 0.
The model includes the effect of convection in the current collector and in the electrode. At steady state the velocity can be derived by putting Eq. 12 into Eq. 5, resulting in At the current collector, the boundary condition is The model is one-dimensional, representing an average behavior of the whole cell area. In the experimental cell the inlet gases enter at the center of the cylindrical cell and flow out in the radial direction. The change in the gas flow rate, V , flowing over a circle of radius r can be written as ) unless CC License in place (see abstract where i app is the applied current density for the electrode , A cell the geometric area of the cell, V m the molar volume of the gas and R cell the radius of the cell. r = R cell /2 is used in this model. The molar fractions at the current collector boundary are calculated based on the inlet molar fractions entering the cell, x i,inlet , using the following relation whereV is the gas flow rate entering the cell. For the case when considering the WGSR, x i,cc are recalculated using a reaction shift coordinate in order to comply with the equilibrium relation (Eq. 13).
At the electrode/matrix interphase, there is no flux of the gas components out from the electrode at steady state, resulting in the following boundary condition.
In the MCEC, the effective conductivity for the solid phase of the Ni electrode is much higher than for the electrolyte phase, indicating that the solid, electron-conducting, phase potential may be set as constant. Thus the potential change is assumed to occur only in the electrolyte phase, which follows the Ohm's law and is given by where i l is the geometric current density in the electrolyte phase, κ eff is the effective electrolyte conductivity, and l is the potential in the electrolyte phase.
The following expression is used to describe the electrode kinetics

RT
[20] Here i 0 0 is the standard exchange current density. The γ coefficients are lumped kinetic parameters that incorporate effects on both reaction order and transfer coefficients for the unknown rate-determining step of the hydrogen evolution reaction. The transfer coefficients α a and α c are unknown and could have been chosen as fitting parameters in the model. However, since most of the measurement data stems from overpotentials below or similar to the thermal voltage (RT/F≈79 mV), we are mostly operating in the linear regime of the exponential terms and the number of data points were deemed not be sufficient for fitting α a and α c . Hence α a and α c are assumed to be 1 in this paper. (Using a linearized kinetics expression with respect to η gave almost identical values for fitting parameters.) In a porous electrode the current density is distributed along the depth of the electrode and the local electrochemical reaction rate is dependent on the local overpotential, η, i.e. the difference from the equilibrium potential. 14 where s is the potential in the solid phase, and E eq is the equilibrium potential for Reaction 1, which is calculated against the reference electrode (i.e. gold wires in equilibrium with the gas mixture of 33/67% O 2 /CO 2 ) and given by where E 0 is the standard equilibrium potential for hydrogen production in the MCEC, and x i,ref is the molar fraction of species i in the reference gas. The value of E 0 is calculated from the literature 15,16 and is listed in Table I. At the current collector/electrode interphase, the boundary condition is At the electrode/matrix interphase, the boundary condition is given by When plotting the simulated polarization curves, we use the following definition for the overpotential of the electrode at the electrode/matrix interphase. [25] The model is solved in Comsol Multiphysics 5.2a (5.2.1.199). A mesh of 101 elements of equal size is used in the model, with the element type set to linear. The SNOPT optimization solver was used to perform a simultaneous parameter fitting of κ eff , i 0 0 , and γ for a set of seven different gas compositions, using a global least square objective function based on the difference between the experimental stationary polarization curves and the corresponding modeled cell polarization.
All parameter values used in this model are presented in Table I.

Results and Discussion
The fitted parameter values for the model, both when including and not including the water gas shift reaction (WGSR), are shown in Table III. The global least square objective function value was about 70% higher when not including the WGSR. It can be seen that the standard exchange current density is similar for the two cases, with the fitted conductivity somewhat higher for the WGSR case. The values of the γ coefficients however differ considerably between the two cases. When the WGSR is not included, the kinetic coefficients of water, carbon dioxide and hydrogen show similar values, indicating that the three reactants have a comparable effect on the rate-determining step of hydrogen evolution reaction. However, when including the WGSR, water shows much stronger influence on the electrode kinetics than carbon dioxide and hydrogen. In the previous study, 10 the reaction mechanisms of hydrogen production in the MCEC, based on the reverse process of hydrogen oxidation proposed by Ang and Sammells 17 and Jewulski and Suski, 18 show kinetic values γ of 0.25 for water, carbon dioxide and hydrogen. The two mechanisms take into account the water-gas shift reaction. However, the fitted kinetic coefficients in this model are not consistent with the theoretical values (0.25) regardless of including the WGSR or not. A possible explanation is that the proposed reaction mechanisms describe the behavior of the planar electrode, which differs from the properties of the porous electrode analyzed in this study. The fitted value of 0.23 for the gas phase porosity when including the WGSR is lower than the value of 0.44 when not including the WGSR. In our model, a higher porosity increases the effective diffusion coefficients, so it seems as if introducing the WGSR in the model reduces mass transport limitations of the cell. Considering the known porosity of 0.54 of the electrode, when not filled with electrolyte, both values are feasible, although for including WGSR the value seems more reasonable in order to get simultaneous percolation of both the gas and electrolyte phases in the electrode. Figures 1-3 show the experimental polarization curves along with the corresponding model values. As can be seen, the visual differences between when including the WGSR and when not including the WGSR are small, and both models give a fair match to the     full range of experimental data. For high current densities, deviations up to 25% are seen between the model and experimental current density values for the same overpotentials. Generally, the model seems to underestimate cell polarization for high current densities, indicating that there are additional polarization losses for high current densities not accounted for in the model. We believe this could be attributed to mass transport effects in the electrolyte phase, which are currently not included in the model. Figures 4-6 show the modeled gas fraction of H 2 , H 2 O and CO 2 at 5000 A · m 2 for a varying inlet molar fraction of the same gas. Regardless of if the WGSR is included in the model or not, all curves exhibit a similar trend, with H 2 increasing in the x-direction (toward the electrode side), and H 2 O and CO 2 decreasing in the x-direction. Only for the case of 10% CO 2 , considerable mass transport limitations seem to be present. This is related to the low diffusion coefficient of carbon dioxide in nitrogen.

Conclusions
In this study a one-dimensional model for the porous Ni electrode of a MCEC based on the Maxwell-Stefan transport equations for mass transfer in the gas phase combined with equations describing the current distribution in the electrolyte phase was developed. The model was fitted to experimental steady state polarization data for various inlet gas compositions of H 2 O, CO 2 and H 2 between 10 and 40%. The model was able to predict cell performance for a broad range of inlet gas compositions. Only in the case of low CO 2 inlet concentrations, mass transport limitations seem occur in the cell at the current levels investigated.
Two separate sets of model parameters were fitted: for when including and not including the WGSR in the model. In both cases it was possible to get a good fit to the experimental data for the same number of fitting parameters, but including the WGSR gave a significantly lower value of the global least squares objective function, and a somewhat more reasonable estimate of the gas phase porosity. It was also seen that the kinetic coefficients get significantly different depending on the choice of model. More experimental data, especially at higher current densities, is needed to shed more light of the role of the WGSR in the MCEC.
The model introduced in this paper provides a better understanding of the porous electrode behavior in the molten carbonate electrolysis cell, and could provide a valuable tool for future studies on, for instance, the behavior of larger electrolysis cells or stack development.

Acknowledgments
The financial support of the China Scholarship Council (CSC) and of StandUp for Energy is appreciated.

List of Symbols
A cell geometric area of cell, m 2 c i concentration of species i, mol · m −3 D i j binary diffusion coefficient between species i and j, m 2 · s −1 D eff i j binary diffusion coefficient between species i and j, m 2 · s −1 E 0 standard equilibrium potential, V E eq equilibrium potential of electrochemical reaction, V F Faraday's constant, A · s · mol −1 i 0 0 standard exchange current density, A · m −2 i app applied current density for electrode, A · m −2 i l geometric current density in electrolyte phase, A · m −2 j loc local current density from electrochemical reaction, A · m −2 K p equilibrium constant for shift reaction n number of electrons N i molar flux of species i with respect to stationary axis, mol · m-−2 · s −1 P pressure, Pa r flow over a circle of radius r = R cell /2, m r shift local reaction rate of shift reaction, mol · m −3 · s −1 R gas constant, J · mol −1 · K −1 R cell radius of cell, m R i reaction rate of species i, mol · m −3 · s −1 S specific surface area, m −1 t time, s T temperature, K u velocity, m · s −1 V gas flow rate entering cell, m 3 · s −1 V m molar volume of gas, m 3 · mol −1 x i molar fraction Greek α symmetry factor in electrochemical reaction γ kinetics coefficients for partial pressure dependencies ε porosity (gas volume fraction) η overpotential, V