Macroscopic multi-dimensional modelling of electrochemically promoted systems

The objective of this work is the construction of macroscopic models for electrochemically promoted catalytic systems, i.e. systems where the catalytic performance is improved by application of potential between the anode and cathode electrodes in the cell. This polarization effect leads to a formation of an effective double layer over the catalytic film due to migration of ‘backspillover’ species from the electrolyte to the working electrode when potential difference is applied in the system. In this paper, we propose a multidimensional, isothermal, dynamic solid oxide single pellet model, which describes the chemical and electrochemical phenomena taking place under polarization conditions. The electrochemically promoted oxidation of CO over Pt/YSZ is used as an illustrative system. The partial differential equation-based 2and 3-dimensional macroscopic models that describe the simultaneous mass and charge transport in the pellet are constructed and solved in COMSOL Multiphysics. The model predicts species coverage on the catalytic surface, electronic and ionic potential curves across the pellet, gas mixture concentration within the reactor and CO2 production rate. Parameter estimation is undertaken so as to provide us with values of parameters, which are necessary for the simulation of the model. Subsequent sensitivity analysis is performed to investigate the effect of the percentage change of each estimated parameter on the CO2 production rate. As it is shown, the reaction rate curves obtained from the current modelling framework are in good agreement with values found in the literature. & 2013 Authors. Published by Elsevier Ltd. All rights reserved.


Introduction
The subject of this work is the development of accurate macroscopic models for electrochemically promoted catalytic systems. Electrochemical promotion of catalysis (EPOC), also referred to as non-Faradaic electrochemical modification of catalytic activity (NEMCA), is the enhancement of catalytic activity and selectivity as a result of an electrochemically controlled migration of "backspilllover" species, i.e. [O δÀ À δþ], from the solid electrolyte to the catalytically active, gas exposed, electrode surface, when current or potential is applied between the two electrodes of the solid electrolyte cell. More specifically, oxygen anions are excorporated from the triple phase boundaries (TPBs), i.e. places where gas phase, catalytic film and electrolyte are in contact) forming BackSpillover Species (BSS) and then the BSS spill over the catalytic surface. Due to this migration, an effective double layer is formed over the catalytic film, which affects the binding strength of the chemisorbed reactants. Fig. 1 illustrates the effective double layer for a system of gaseous reactants (CO,O 2 )/Metal (working electrode)/YSZ (solid electrolyte). The EPOC phenomenon was first observed in the early 1980s by Stoukides and Vayenas (1981) and has since been extensively studied experimentally and theoretically (Ladas et al., 1991;Metcalfe, 2001;Poulidi et al., 2011; Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/ces et al., 1990;Xia et al., 2010). However, few modelling studies have been performed describing the addressed effect, which were primarily focused on the catalytic surface kinetics. To the best of our knowledge, a systematic multi-dimensional modelling framework, which simulates the chemical and electrochemical phenomena taking place in electrochemically catalyzed systems and will allow robust system design, has not been formulated yet. Stoukides and Vayenas (1981) have observed that the nature of a catalyst deposited on a solid electrolyte can be altered in a dramatic manner, when voltages are externally applied to the solid electrolyte cell. The study included ethylene oxidation on polycrystalline silver as catalytic reaction, at atmospheric pressure and temperatures close to 400 1C. It was observed that O 2 À pumping from the electrolyte to the catalytic surface was taking place by applying voltage between the cathodic and anodic electrodes. This O 2 À pumping phenomenon was found to result in an increase in the rate and selectivity of the heterogeneous catalytic reaction in a considerable and reversible manner. By the late 1980s, NEMCA was demonstrated for several catalytic reactions and was not limited to a specific solid electrolyte or catalyst (Vayenas et al., , 1989. It was shown that the work function of a catalyst can be altered when oxygen anions migrate from the solid electrolyte to the catalytic surface resulting in significant changes in the catalytic activity and as much as 4000% increase in the surface reaction rate. An enhancement factor, Λ (Λ ¼ Δr=ðI=2FÞ), was introduced to report whether the behaviour of the system was Faradaic or not . Here, Δr describes the polarization-induced change in the catalytic rate and I/2 F is the rate of supply/removal of oxygen anions (Faradaic rate) when electrical current, I, or a potential difference between the working and the reference electrodes is applied. A value of unity for Λ means that the behaviour of the system is purely Faradaic. Greater values signify non-Faradaic behaviour. Vayenas et al. (1990) have shown that the catalytic activity is exponentially dependent on the catalyst work function, which is in turn affected by polarization conditions in the cell.
The effect of EPOC was further studied giving insights of the surface dynamics of a catalytically electropromoted system (Ladas et al., 1991;Vayenas and Bebelis, 1997;Yentekakis et al., 1994). It was, moreover, found that the NEMCA effect was not limited to pure catalytic films, but it was also observed when dispersed (Marwood and Vayenas, 1998) or metal oxide (Wodiunig and Comninellis, 1999) catalysts were utilized. Ladas et al. (1991) established that just a small coverage of ions (o5% of a monolayer) can change the work function of the catalyst without affecting the coverage of the rest chemisorbed species, while Yentekakis et al. (1994) observed that for Na coverages lower than 0.04, sodium was acting as promoter by increasing the surface reaction rate up to 600% and for greater than 0.04 Na coverage values, sodium acted as poison decreasing the rate by up to 90% due to the formation of a surface complex between CO, Na and Pt. Vayenas and Bebelis (1997) using X-ray photoelectron spectroscopy provided significant evidence that the normally chemisorbed oxygen, i.e. oxygen in the gas phase which is dissociatively adsorbed on the catalytic surface, was greatly more reactive than the backspillover one. Marwood and Vayenas (1998) studied NEMCA in the case of a dispersed catalyst, Pt, on Au. It was observed that the NEMCA had the same features as in the earlier studied pure metal catalysts. Furthermore, Wodiunig and Comninellis (1999) used a metal oxide catalytic film, RuO 2 , and found that the O 2 desorption activation energy was getting lower by increasing the potential difference in the cell. This behaviour was completely reversible under current interruption.
Recently, Pliangos et al. (2000), Vayenas et al. (2003) as well as other research groups (Poulidi et al., 2007(Poulidi et al., , 2011Xia et al., 2010) have investigated the electrochemical promotion phenomenon, using novel configurations and bringing new challenges in the field of modern electrochemistry. Pliangos et al. (2000) investigated the NO reduction via dry impregnation with NaOH and found that the reduction rate was further increased by up to a factor of 20 under polarization conditions. Vayenas et al. (2003) have shown using scanning tunneling microscopy that the BSS can diffuse on the catalytic surface for greater than mm distances and thus the NEMCA effect could take place in commercial nm sized catalytic films. Poulidi et al. (2007) used a novel wireless configuration in a dual chamber reactor. More specifically, they have used a mixed protonic-electronic conductor as electrolyte to supply hydrogen promoting species on the catalytic surface (reaction side chamber) by manipulating the hydrogen chemical potential difference throughout the electrolyte. The adjustment of the chemical potential of hydrogen across the electrolyte was performed by utilizing hydrogen, helium and oxygen in helium as sweep gas (sweep side chamber). It was observed that the reaction rate could be increased by up to a factor of 1.6 only in the case of hydrogen as sweep gas and under certain conditions. This enhancement effect was completely reversible. Later, Poulidi et al. (2011) continued their wireless configuration studies utilizing a mixed oxygen ionic-electronic conductor as electrolyte in several experimental configurations. In this case, it was found that although there was an enhancement of the catalytic activity under oxygen utilization as sweep gas, the reaction rate did not return to its unpromoted state. Another novel configuration was developed by Xia et al. (2010). In this study, CO combustion was performed on a bipolar configuration of nano-dispersed Pt particles deposited on YSZ. Application of current was found to increase the reaction rate as much as 500 times. This non-Faradaic promotion effect was irreversible after current interruption.
As mentioned earlier, despite the great potential of modelling EPOC, as it will allow robust system design and thus the establishment of industrial-level applications, very few modelling studies have been performed so far. Vayenas and Pitselis (2001) developed a steady-state, 1D surface reaction-diffusion model to simulate NEMCA for a porous catalyst deposited on a solid electrolyte. In that study, BSS was the only species taken into account and the TPB was considered to represent the interface layer (with negligible thickness) between the porous catalytic film and the solid electrolyte. Furthermore, the BSS could be excorporated from the electrolyte and either diffuse through the catalytic film or react with the gas phase species. This model was used in Presvytes and Vayenas (2007) to study a system consisting of a semi-spherical catalytic particle interfaced with a larger YSZ one. The TPB, and thus the reaction zone, was extended over the entire thickness of the catalytic particle. An effectiveness factor was used to determine the catalytic performance. It was found that the catalytic particle geometry has a minor effect on the catalytic performance. Foti et al. (2002) formulated a model to describe the reaction rate transients of an electrocatalytically promoted system. Fast diffusion of BSS, which could form on the catalytic surface if and only if there were free active surface sites, on the catalytic surface as well as a first order BSS consumption rate were considered. It was found that the BSS coverage on the catalytic surface was  (Vayenas and Bebelis, 1999). exponentially dependent on the time of polarization and that the steady state value of the BSS coverage, as well as its accumulation on the catalyst rate were increasing functions of the current application. Leiva et al. (2008) used Monte Carlo simulations in the Grand Canonical ensemble to study NEMCA. They analyzed the effective double layer to explain the modification of the catalytic surface due to BSS migration. It was found that Na þ and O δ À promoting species had different effect on the lattice of the catalyst. The former seemed to generate an expansion of the catalytic lattice as a result of transfer of electrons to the metal, while the latter was found to induce a compression of the lattice due to the opposite effect. Recently, we have performed a modelling study (Bonis et al., 2011) to describe electrochemical promotion. We developed a 2-D macroscopic model to simulate reaction-diffusion kinetics as well as electrostatic and electrochemical phenomena for the electrocatalytically promoted CO oxidation. However, only qualitatively results were obtained due to the lack of reliable values of the parameters used in that study.
The main scope of this work is the development of a systematic validated macroscopic modeling framework of the chemical end electrochemical processes taking place in a solid oxide single pellet under the effect of electrochemical promotion. The electrochemically promoted oxidation of CO over Pt/YSZ is used as an illustrative system. Such a model can be a valuable tool for the design of "industrial" scaled up systems such as exhaust gas treatment, hydrogenation of organic compounds and product selectivity modification and solid oxide fuel cells. To accomplish this, the present study focuses on the formulation of a multidimensional, isothermal, dynamic solid oxide single pellet model by coupling mass transport (surface reaction-diffusion) with charge transport and electrochemical phenomena (electrochemical reactions taking place at the TPBs of anode and cathode electrodes). Parameter estimation is performed using experimental data from the literature and a validated model is obtained. Subsequent sensitivity analysis of the estimated parameters is performed to assess the effect of the percentage change of each parameter on the reaction rates. Finally, the effects of the partial pressures of the gaseous reactants as well as of the temperature on the reaction rate are further investigated and conclusions from the above analysis are discussed.

The chemical process
The electrochemically promoted oxidation of CO over Pt/YSZ will be used to demonstrate the detailed model constructed in this work, since CO oxidation on platinum is one of the most extensively studied heterogeneous catalytic systems. Fig. 2 illustrates the surface dynamics of CO oxidation.

Reaction mechanism
Several studies have shown that the CO oxidation on Pt follows the Langmuir-Hinshelwood mechanism (Kaul and Wolf, 1985;Palmer and Smith, 1974). This mechanism involves a surface reaction between adsorbed CO and oxygen atoms forming CO 2 (Fig. 2). Once CO 2 is produced, it desorbs immediately and its readsorption on the Pt surface is assumed to be negligible (Matsushima, 1978). Chemisorption of O 2 and of CO on Pt occurs through dissociation (Ducros and Merrill, 1976) and without dissociation (Heyne and Tompkins, 1966), respectively. The corresponding mechanism is described below (reaction (1)-(3)) (Kaul et al., 1987): When potential is applied between the anode and the cathode of the solid oxide single pellet, backspillover species are formed at the anodic TPBs and "migrate" to the Pt surface (Fig. 3). These species can either react with the adsorbed CO forming CO 2 (reaction (4)) or desorb in the gas phase as O 2 (reaction 5) (Vayenas and Pitselis, 2001).

Surface mass balances
The surface mass balances for atomic O 2 , CO and BSS are described by Eqs. (6)-(8), respectively. It should be noted that the surface reaction is not assumed to be the rate determining step and that the adsorption/desorption reactions are not considered to be at equilibrium. The conservation of active catalytic sites, θ S , can be expressed by Eq. (9). Hence, where R i is the overall reaction rate for species i, θ j the coverage of species j on the catalytic surface, r i the reaction rate for reaction i (reactions (1)-(5)), and r elec is the Faradaic production rate which is described below (see Eq. (40)). The individual reaction rates are defined as follows (Kaul et al., 1987): Fig. 2. Schematic presentation of CO oxidation on platinum surface. Fig. 3. Schematic presentation of BSS formation and "migration" on Pt.
where k i is the rate constant for reaction i, C A O 2 and C A CO the gaseous concentration (the superscript A denotes the anode) of O 2 and CO, respectively, which can be calculated by the ideal gas equation (Poling et al., 2001): where P j is the partial pressure for species j in the gas phase, R the ideal gas constant and T the temperature in the gas phase. The terms (1 À θ O ) 2 and (1 À θ BSS ) 2 are incorporated in the denominators of oxygen's and backspillover species' desorption rates in order to limit the saturation coverage (Gorte and Schmidt, 1978;Kaul et al., 1987). The rate constants for adsorption for both oxygen and carbon monoxide are assumed to be sticking coefficient-dependent, while the ones for desorption and for the surface reaction are considered to follow an Arhenius expression (Kaul et al., 1987). Thus, the individual rate constants can be given from: where S j is the sticking coefficient, M j the molecular weight of species j, T s the temperature of the catalytic surface and N s the concentration of active sites on the platinum surface. Furthermore, where k o,i and E A,i are the pre-exponential factor and the activation energy for each reaction rate, respectively. The sum of the consumption rate constants of backspillover species can be taken to be equal to 10 À 2 s À 1 (Vayenas and Pitselis, 2001). Thus, 3. The electrochemical process

Electrochemical reactions
At the TPBs of the cathode, the electrochemical reduction of O 2 takes place: At the anodic TPBs, three electrochemical reactions can occur, i.e. oxygen anions can either be excorporated from YSZ as O 2 and/or BSS, or react with CO giving CO 2 . Hence, The rate of the cathodic reduction of O 2 (reaction (17)) is equal to the sum of the rates of the anodic electrochemical reactions ( (18)-(20)), so that the number of electrons consumed by the electrochemical reactions at the cathode TPBs is equal to the ones produced from the electrochemical reactions at the anode TPBs.

Equilibrium potentials
For the cathodic electrochemical reaction (17), the electrochemical potentials of the reacting species have to be balanced at equilibrium (Bard and Faulkner, 2001;Chen et al., 2010). Hence, where μ C O 2 and μ C O 2 À YSZ are the chemical potentials (the superscript C denotes the cathode) of O 2 and O 2 À YSZ , respectively, F is the Faraday constant, Φ C el and Φ io are the local equilibrium potentials of the cathodic electrode and the electrolyte phase, respectively. The equilibrium potential, E eq , is defined as the equilibrium electric potential difference between two phases. Thus, rearranging Eq. (21), the equilibrium potential can be given by: Similarly, the equilibrium potentials E eq,A1 , E eq,A2 and E eq,A3 related to the anodic electrochemical reactions (18)-(20), take the form: where Φ A el is the local equilibrium potential of the anodic electrode and m i is the chemical potential of species i. The superscript A denotes the anode. The factor of 2 in the denominator represents the number of electrons, which are produced and/or consumed by each electrochemical reaction.

Thermodynamic open circuit potential
The thermodynamic open circuit potential, V OC , also known as Nernst potential at open circuit conditions, is expressed in Eq. (26) in an analogous fashion to the one presented in Ho et al. (2009).
Combining Eq. (26) with Eqs. (22)-(25) and assuming that the gas phase mixture behaves ideally (μ i ¼ μ o i þ RT ln P i ) and that the chemical potentials of oxygen anions at the two sides of the electrolyte are equal (Ho et al., 2009) ,we derive the following expression for V OC : where μ A BSS is the chemical potential of the BSS, and V o OC the ideal Nernst potential (under standard conditions) which is expressed as: where the superscript o denotes standard conditions.

Current density distribution
The charge transfer between the electronic and ionic current, due to the electrochemical reactions at the anodic and cathodic TPBs, can be computed by the Butler-Volmer equation. Thus, the current density distribution due to the electrochemical reduction of oxygen (reaction (17)) taking place at the cathodic TPB is expressed as (Suzuki et al., 2008;Wang et al., 2007): where J C 0 is the exchange current density of the cathode, α C the cathodic charge transfer coefficient, n e the number of electrons transferred in the cathodic electrochemical reaction and η C the overpotential of the cathode.
The cathodic reduction of O 2 is taking place three times, thus it is considered as three different reactions and the parallel electrical circuit analogy is applied (Achenbach, 1994;Andreassi et al., 2009), hence the factor of 3 in the above expression.
The exchange current density follows an Arrhenius expression and for the cathodic electrochemical reaction (reaction (17)) it is given by (Andreassi et al., 2007): where γ C is the pre-exponential coefficient, P ref the reference pressure in the reactor, i.e. the total gas pressure in the reactor, and E C A the activation energy of the cathode. The overpotential of the cathode is defined as (Bove and Ubertini, 2006;Tseronis et al., 2008Tseronis et al., , 2012: At the anodic TPB, three electrochemical reactions take place. Hence, the total current density distribution of the anode can be expressed using the parallel electrical circuit analogy as (Achenbach, 1994;Andreassi et al., 2009): where J A 1 , J A 2 and J A 3 are the current densities due to the electrochemical reactions (18)-(20), respectively, given by the following expression: where J A 0;i is the exchange current density of the anode for each electrochemical reaction i, α A the anodic charge transfer coefficient and η A the overpotential of the anode.
The exchange current density of the electrochemical reaction (19) can be described by a similar correlation to the one used for the electrochemical H 2 oxidation and has the following form (Costamagna et al., 2004;Tseronis et al., 2008Tseronis et al., , 2012: where γ A,2 is the pre-exponential coefficient of electrochemical reaction (2) (reaction (19)), and E A A the activation energy of the anode.
The expression for the exchange current density for reaction (18) follows: The power of À 0.25, comes from the stoichiometric analogy between reactions (17) and (18). The exchange current density for reaction (20) is expressed for first time in this study and is assumed to be partial pressure independent due to the unknown nature of BSS. Hence, The overpotential of the anode is defined as (Bove and Ubertini, 2006;Tseronis et al., 2008Tseronis et al., , 2012: where Φ A el and Φ io are the local electronic and ionic potentials at the anode side, respectively.

Model of electrochemically promoted CO oxidation
The well mixed (CSTR) reactor used in this study is of "singlepellet" type and the 3-D computational domain of the single pellet are depicted in Fig. 4. The temperature in the reactor is considered constant. The pellet consists of the electrolyte (YSZ) and three electrodes, the anode-working electrode (Pt), and the cathode counter and reference electrodes (Au), with surface areas of 2 cm 2 , 1.2 cm 2 and 0.1 cm 2 , respectively. The thickness of the electrolyte is 400 μm. The electrolyte and the electrodes are all exposed to the reacting gas mixture. The gas mixture at the inlet of the reactor consists of O 2 and CO while helium is used as carrier gas.
This solid oxide single pellet model can be used to simulate the mass and charge transport as well as the chemical and electrochemical phenomena taking place on the catalytic surface and the triple phase boundaries respectively. COMSOL Multiphysics 3.5a (COMSOL, 2008) is a commercial CFD software (FEM) which is used here to construct and simulate the proposed model. The model computes species coverage profiles at the catalytic surface as well as electronic and ionic potential throughout the pellet, gas mixture concentration in the reactor as well as CO 2 production rate.
Charge transfer occurs in the solid oxide single pellet due to the applied potential difference between the anode and the cathode. Electrochemical reactions take place at the TPBs of the anode and the cathode, leading to a transfer of electronic to ionic current and the vice versa. BSS is formed at the TPBs of the anode and spills over the catalytic surface where CO oxidation is taking place. A 2-D version of the model is presented here. The 2-D computational domain that is used for the solid oxide single pellet is depicted in Fig. 5. It represents an X-Z cross section of the 3-D single pellet presented in Fig. 4. The Pt and Au electrodes are considered as boundaries on the electrolyte rather than finite regions.

Model assumptions
The main assumptions that have been made for this 2-D macroscopic model are listed below: The thickness of the anodic and cathodic electrodes is considered to be negligible since it is too small compared with the thickness of the electrolyte. Thus, the electrodes are 0dimensional boundaries of the electrolyte domain.
The TPBs of the anode and the cathode are respectively represented by the perimeters of the anodic and cathodic electrodes.
The gaseous species around the single pellet are considered to be well-mixed. Moreover, the gas mixture is assumed to behave as an ideal gas.
The overall pressure of the gas phase is considered constant throughout the reactor.
The temperature on the catalytic surface as well as in the solid oxide single chamber and throughout the reactor is assumed constant.
The operating potential difference throughout the pellet is considered to be constant, assuming that both electrolyte and electrodes are good charge conductors.
The electrochemical reactions are considered to take place only at the triple phase boundaries.
The model accounts only for catalytic reactions on Pt while the catalytic activity on Au (in the presence of Pt), can be neglected at as high temperatures as 600 K.
Only the diffusion of BSS is taken into account for the mass transport over the Pt surface, since for the other species adsorption and desorption will dominate, and the relevant coefficient is considered to be temperature and pressure independent.
The chemical potential of BSS is assumed to be equal with the one of the O 2 À on the surface of YSZ (μ BSS ¼ μ O 2 À YSZ ). The value of the "standard" chemical potential μ O 2 À YSZ was found to be equal with À 236.4 kJ mol À 1 (Bessler et al., 2007;Vogler et al., 2009).

Mass transport equations on the catalytic surface
On the non-porous metal phase, the mass conservation species i, has the form (Vayenas and Pitselis, 2001): Eq. (38) can be written in terms of the coverage, θ i , of species i, as: where R′ i is the overall reaction rate of species i (expressed in 1/s units) which is described by Eqs. (6)-(8). The Faradaic production rate, r elec , is described as: where J A 3 describes the current density of anode due to BSS formation (reaction 20).

Charge transport equations in the pellet
The charge conservation in a non-porous medium j (j¼io, el (electronic, ionic)), can be described by the following equation (Cheddie and Munroe, 2006;Tseronis et al., 2012): where ρ j is the charge density, Q j the charge source term and J j the current density. The current density can be described by Ohm's law as: where s j is the charge conductivity and Φ j the potential. Combining Eq. (41) with Eq. (42) yields: Electronic and ionic charge transfer occurs in the solid oxide single pellet due to polarization. The governing equations describing the charge transport at the electrolyte and the electrode domains, considering no charge source, are described as follows: where ρ io and ρ the ionic and electronic potentials respectively. The superscripts A/C denote the anode/cathode.

Gas phase modelling
As mentioned above, the gas phase around the single pellet is considered to be well-mixed and to behave as an ideal gas. The partial pressures of the gas phase in the reactor are described by Eqs. (46)-(49): where P i is the partial pressure of species i, F d the volumetric flowrate of the gas mixture in the inlet/outlet of the reactor and R i the overall gaseous production (i ¼CO 2 ) and consumption (i¼CO, O 2 ) rate due to both chemical and electrochemical reactions. Carbon dioxide is produced by the catalytic surface reactions ( (3) and (4)) and by the electrochemical reaction (19). Carbon monoxide is produced/consumed by the catalytic reaction (2) and electrochemical reaction (19). Oxygen is produced/consumed by the catalytic reactions ( (1) and (5)), by the electrochemical excorporation of oxygen anions at the anodic TPBs (18) and by the electrochemical oxygen reduction (17). These phenomena are summarized in the expressions: where r i is the reaction rate of reaction i as in Eqs. (10)- (12), N S the concentration of the active surface sites on the catalytic surface, W the width for both the anode and cathode, J i /n e F the Faradaic rate resulting by the electrochemical reactions i, L the length of the anodic electrode and L count the length of the counter electrode in the cathode.

Boundary conditions
The boundaries (B) of the 2D computational domain are shown in Fig. 6. There are in, total 8 boundaries and 8 points, 4 of which, P1, P3, P6 and P8, are considered as TPBs. It should be noted that the working electrode in the anode (Pt) is represented by boundaries B1 and B2 in this computational domain while the counter electrode in the cathode (Au) is represented by B6 and B7.

Mass transport.
At the edges of the anodic electrode (points P1 and P3 in Fig. 6) the no flux condition is imposed for all species. Hence, P1 and P3 : n Â ðÀD i ∇θ i Þ ¼ 0; where n is the unit normal vector.

Charge transport.
The circuit is closed at points P7 and P2 (see Fig. 6). The electronic potential is fixed to the value of the operating potential Φ pellet at P7 and is equal to 0 at P2. At the points P6 and P8, the cathodic reduction of O 2 takes place and leads to a transfer of electronic to ionic current. Thus, P6 and P8 : where J c is computed from the modified form of the Butler-Volmer (Eq. (29)). At the points P1 and P3, transfer of ionic to electronic current occurs, due to the electrochemical reactions. Hence, P1 and P3 : here J A is computed by Eqs. (32) and (33). At all the remaining boundaries and points, no flux conditions are imposed.

3-D macroscopic model
The 3D computational domain of the solid oxide single pellet is illustrated in Fig. 4. Similarly to the 2D model, the electrodes of the anode and the cathode are considered as boundaries of the electrolyte rather than finite regions. The mass and charge transport phenomena are, here too, described by Eqs. (39), (44) and (45). The gaseous consumption/production rates are given by the surface integrals: 4.3.1.2. Charge transport. The electronic potential is fixed to the value of the operating potential Φ pellet at point P6 and to 0 at P5.
Transfer of electronic to ionic current occurs at E6 À7 and E13 À 14. Hence: E6 À7 and E13 À 14 : Ionic charge is transferred to electronic at E1 À 4: At all the remaining boundaries and edges the no flux condition is imposed for both ionic and electronic charge transfer.

Numerical solution approach
The finite elements method (FEM) is used for the solution of the set of partial differential Eqs. (39,44,45) including (53)-(59) implemented through COMSOL Multiphysics 3.5a. The 2-and 3-D computational domains were discretized in 6568 triangular and 29,438 tetrahydral triangular elements, respectively which were deemed sufficient to produce mesh-independent solutions. The resulting sets of nonlinear algebraic equations were solved using the Direct (UMFPACK) solver in COMSOL.
The charge transport in both the anode and cathode electrodes as well as the mass balances on the catalytic surface were implemented by using the weak form boundary PDEs. Moreover, the expressions for each species partial pressure in the reactor (Eqs. (46)-(48)) were implemented by utilizing the COMSOL facility of "integrated coupling variables". In Table 1, the chosen model parameters are presented. The electronic charge conductivities are considered to be temperature dependent and are given by Eqs. (68) and (69)  The ionic charge conductivity is given by (Ferguson et al., 1996): Table 1 Solid oxide single pellet model parameters.

Results and discussion
The system operating conditions used in the model are tabulated in Table 2. The pellet is considered to be preheated at temperature T ¼645.15 K and the operating potential Φ pellet is assumed to be constant throughout. Extensive parameter estimation studies were performed, using the 2-D model, for parameters, with unknown values, which are of crucial importance for modelling EPOC systems. Furthermore, both the 3-and the 2-D model were validated against experimental data from the literature (Yentekakis et al., 1994). Furthermore, the 2-D model was exploited for sensitivity analysis as well as parametric investigation studies.

Parameter estimation
The reaction kinetic constants for this system are not available in the literature. Although the catalytic system for the oxidation of CO on Pt has been extensively studied (Campbell et al., 1980(Campbell et al., , 1981a(Campbell et al., , 1981b, a direct extension of catalytic kinetic data to the electrochemically promoted system is not plausible: polarization changes both the, steady state and the dynamic, system behaviour. This observation has been the backbone of all NEMCA studies to date. The underlying cause for this change is the electrochemically induced alteration of all reaction and transport phenomena taking place stemming from (i) the addition of reactions (4) and (5) to the Langmuir-Hinshelwood mechanism for the CO oxidation (reaction (1)-(3)) and (ii) the rate change for reactions ((1)- (3)). In fact, the former is less significant, as for this system it has been reported that the rates of reactions ((4)-(5)) are relatively low and constrained by Eq. (16). Thus, the kinetic parameters of the closedcircuit system are expected to differ significantly from the ones of the open-circuit system. This is especially so since the reported bounds for the given open-circuit kinetic values are quite broad (cf. Kaul et al., 1987).
For this reason, we performed a parameter estimation study, using our 2-D model in conjunction with reported closed-circuit experimental data of overall CO 2 reaction rates for a range of CO partial pressures available in literature (Yentekakis et al., 1994). The operating potential was 700 mV, the partial pressure of the oxygen at the reactor inlet was 5.8 kPa, and the temperature of the gas mixture in the reactor was 372 1C.
The geometrical aspects (as in Table 3) of our model used for the estimation studies are considered to be typical for this kind of experiments (Vayenas et al., 1989;Yentekakis and Bebelis, 1992;Yentekakis et al., , 1994. The system electrodes were reported to have thickness of 5 μm, which we considered to be negligible compared with the 0.4 mm thickness of the electrolyte (YSZ), hence they were modeled as "flat" surfaces (see Figs. 4 and 5).
The kinetic parameters that were estimated through our model were: the sticking coefficient of O 2 (S O2 ) and CO (S CO ) on Pt (reactions (1)-(2), Eq. (14)), O 2 and CO desorption activation energies (E A, À 1 and E A, À 2 respectively, reactions (1)-(2), Eq. (15)), the CO oxidation activation energy (E A,3 , reaction (3), Eq. (15)), the BSS surface reaction rate constants (k 4 or k 5 , reactions (1)-(2), Eq. (15)) as well as the pre-exponential coefficients for the three electrochemical reactions of the anode (γ A,1 , γ A,2 and γ A,3 , reactions (18)-(20), Eqs. (34)-(36)) utilized in Butler-Volmer expressions (see Table 1). The parameter estimation problem is set up as a nonlinear least squares problem: where n wxp is the number of experimental points,R CO 2 the mean overall rate for CO 2 production, i.e. the average rate due to chemical and electrochemical processes, which is a function of the inlet pressure of CO,R CO 2 ðP in CO Þ, as defined in Eq. (50), and R CO 2 the experimental CO 2 production rate for the same P in CO . The model Eqs. (16), (39), (44)-(59) provide the equality constraints for this optimization problem.
We assume that the NEMCA effect affects only the activation energies of the catalytic reaction rates (which is reasonable since closing the circuit will in general increase the catalytic activity by affecting the activation barriers), while the values of the reaction rate constants are taken from equivalent open-circuit experiments (Kaul et al., 1987). On the other hand, the activation energies for the anodic electrochemical reactions are assumed to be equal to a value taken from literature (see Table 1) and here we have selected to estimate the anode pre-exponential coefficients due to the unknown nature of BSS and the assumption made for the exchange current density of reaction (20) (see Eq. (36)). Finally, the activation energy of the cathodic electrochemical reaction is taken from literature (Chaisantikulwat et al., 2008) (see Table 1) as the oxygen reduction taking place at the cathodic TPBs is typical for such systems. The parameter estimation task was performed with the use of the 2-D model (whose computational domain is illustrated in Fig. 5) for computational efficiency. The computed parameters were subsequently introduced in the 3-D model for validation purposes.
It is worthwhile to note here that the limited number of experimental data to only one polarization potential, do not allow to investigate the parametric effect of voltage and also result in identifiability issues, noisy objective functions and multiple extrema.
Parameter estimation studies for electrochemical systems in the context of SOFCs have been undertaken with the use of various optimisation methodologies (e.g. Shearing et al., 2010).
In this work to overcome the problem of convergence to local minima, we have solved the estimation problem using Simulated Annealing (SA) (Kirkpatrick et al., 1983), a stochastic optimization method which, with appropriate restart techniques, can help to avoid convergence, merely, to the nearest local optimum point. On the contrary, this algorithm is expected to probabilistically converge to the neighborhood of the global optimum. A refining step using a gradient-based, deterministic strategy, Sequential Quadratic Programming (SQP) (Gill et al., 1981) implemented through the fmincon function in MATLAB, is subsequently performed using as initial guess the result from SA. In all cases fmincon converged Table 2 Model operating conditions.

Domain Length (cm) Width (cm) Other characteristics
Anode 1 2 Area: 2 cm 2 Electrolyte 1 2 Volume: 0.0008 cm 3 Cathode (counter) 0.6 2 Area: 1.2 cm 2 monotonically within 7-9 iterations. It is possible, although improbable, that the attractor of the global minimum is weak and SQP will converge to a local minimum of higher objective function value. Here the best of the SA/SQP solutions is taken as the result of our estimation problem. The values of the parameters identified using the procedure outlined above are delineated in Table 4. As it can be seen all sticking coefficients and activation energies are within the upper and lower limits reported in (Kaul et al., 1987). Furthermore the computed CO desorption activation energy and the O 2 sticking coefficient are lower than the reported ones for the catalytic system (Kaul et al., 1987) which follows the NEMCA theory suggesting weaker bonds of the adsorbed species due to the polarization effect .
It should be noted here that the kinetic parameters computed though the fitting procedure, although they do illustrate the BSS effect, since they deviate from the reported open circuit values (Kaul et al., 1987), they are by design constant hence independent of the BSS coverage variation. Nevertheless, it has been reported in (Brosda and Vayenas, 2002) that the kinetic constants are affected by the coverage variations of the BSS species and thus by the catalyst work function. However, the available experimental data, are limited to one temperature and one operating potential, hence this effect is not currently included in our model.
In Fig. 8, the good agreement between model predictions using the computed parameters and experimental data is depicted. Moreover, in Fig. 9, the predictions of the 3D model, using the estimated parameters in Table 4 compared with the experimental data are illustrated also showing an excellent agreement. This signifies that the current system can be sufficiently (and efficiently) represented by the 2-D model.
Uncertainty in the parameter estimation is investigated through sensitivity analysis.

Sensitivity analysis
In order to quantify the effect of variations in the estimated chemical and electrochemical kinetic constants, a sensitivity analysis was carried out, using the 2-D model, by introducing a percentage change ( À50 to þ100%) in the values of these parameters. As it is shown in Table 5, changes in the sticking coefficients and activation energies have a strong effect on the CO 2 production rate, while the rate constants of the reactions involving BSS as well as the exchange current density pre-exponential coefficients have almost no effect on the CO 2 production rate as expected since BSS mainly forms the effective double layer which affects the binding strength of the other adsorbates, while the exchange current density is related to the Faradaic contribution (as opposed to NEMCA) which is in great agreement with research studies suggesting that the enhancement of CO 2 production rate comes mainly from the modification of the catalytic surface leading to an increase in the surface reaction rate (NEMCA effect). This is corroborated in Fig. 10 which illustrates the effect of P in CO on the enhancement factor Λ. Values of Λ of the order of 10 3 (for P in CO 41) are in good agreement with the ones found in literature (Vayenas and Bebelis, 1999) for this system and show how strong the non-Faradaic effect is compared with the Faradaic one. It should be noted that the increase of O 2 desorption Table 4 Values of estimated parameters.

Parametric study
We have also performed parametric studies, using our 2-D model, in order to investigate the effect of operating conditions on CO 2 production rate.
Effect of P in CO . Figs. 16 and 17 depict the effect of P in CO on the CO 2 production rate due to non-Faradaic and Faradaic effects respectively, at constant values of Φ pellet ¼700 mV and P in O 2 ¼ 5:8 kPa. The rate curve presented in Fig. 16 exhibits a behaviour characteristic for Langmuir-Hishelwood kinetics. A "volcano-type" behaviour is observed in the CO 2 production rate at relatively high P in O 2 values due to the NEMCA effect which is in good agreement with literature (Yentekakis et al., 1994), while an "S-type" behaviour is favoured in the CO 2 production rate due to the Faradaic effect (Fig. 17). We can also observe a sharp increase for both Faradaic and non-Faradaic rates at low P in CO , i.e. values between 1 and 2 kPa. Effect of P in O 2 . The effect of P in O 2 on non-Faradaic and Faradaic reaction rates at constant values of P in CO ¼ 1.5 kPa and Φ pellet ¼700 mV is illustrated in Figs. 18 and 19 respectively. In this case, "S-type" behaviour is favoured in the CO 2 production rate due to both NEMCA and Faradaic effects. In Fig. 18 a CO 2 non-Faradaic rate increase is observed for values of P in O 2 in the range of 5-10 kPa, while the rate seems to be reaching a plateau for relatively high P in O 2 . On the other hand, a relatively small rate increase due to the Faradaic effect is observed even for higher values of P in O 2 . Effect of temperature. In the case of constant P in O 2 ¼5.8 kPa and varying P in CO , increasing the temperature leads to a minor increase of the electropromoted rate for values of P in CO lower than 1.5 kPa, Fig. 10. Effect of inlet partial pressure of CO on enhancement factor Λ. Fig. 11. Effect of inlet partial pressure of CO on CO 2 production rate for various values of O 2 desorption activation energy, E À 1 . Fig. 12. Effect of inlet partial pressure of CO on CO 2 production rate for various values of CO desorption activation energy, E À 2 . while the opposite effect is observed for higher P in CO . Likewise, in the case of constant P in CO ¼1.5 kPa and variable P in O 2 , an increase in temperature leads to an increase of the CO 2 non-Faradaic production rate for P in O 2 /P in CO o6 (excess of CO), while a plateau is reached for P in O 2 /P in CO 46. Increasing the temperature in all cases leads to an increase of the CO 2 production rate due to the Faradaic effect.

Conclusions
The objective of this work was to construct a multi-dimensional, isothermal, dynamic solid oxide single pellet macroscopic model to describe the chemical and electrochemical processes that occur in a solid oxide pellet under the effect of electrochemical promotion. The proposed model comprises PDEs for mass and Fig. 15. Effect of inlet partial pressure of CO on CO 2 production rate for various values of surface reaction between O 2 and CO activation energy, E 3 . Fig. 16. Effect of inlet partial pressure of CO on CO 2 production rate under NEMCA effect for various temperature values.    charge balances describing the chemical and electrochemical phenomena taking place at the catalytic surface and the TPBs of the anode and cathode, and couples the operating voltage to adsorbed species coverage and surface reaction rates. The simulation of the model was performed utilizing COMSOL Multiphysics, a commercial CFD software which employs the Finite Elements method to solve the resulting set of PDEs simultaneously. Additionally, the facility of "integrated coupling variables" incorporated in this software was used for the implementation of each gas species partial pressure expressions. The model allows the prediction of species coverage on the catalytic surface, electronic and ionic potential profiles throughout the pellet, gas mixture concentration within the reactor and CO 2 production rates.
Since parameters such as sticking coefficients and activation energies for the electropromoted catalytic reactions are expected to be different from the open-circuit ones and electrochemical parameters such as the pre-exponential exchange current density coefficients are largely unavailable, we undertook a parameter estimation study. To accomplish this, a combination of stochastic and deterministic optimization studies were undertaken utilising a tailored 2D version of our model (through the use of the scripting facility of COMSOL Multiphysics) to match model predictions with experimental data available in the literature (Yentekakis et al., 1994). Possible identifiability issues due to the limited number of available data were addressed through multiple restarts for the stochastic optimiser. The estimated parameters were within literature reported upper and lower bounds. The resulting 2-and 3-D models predict physically realistic trends.
It should be noted here that an increase in P in CO incurs a decrease in the BSS coverage as expected from the exchange current density expressions. However, the effect of P in O 2 is not explicitly quantifiable from the available experimental data, hence further experimental information would be needed to quantify the effect of BSS coverage changes on the catalyst work functions and hence on the kinetic constants.
Sensitivity analysis for the estimated parameters was performed in order to obtain insights on the reliability of the parameter estimation problem and to quantify the effect of each parameter on the reaction rate. As a result, we observed that changes in the chemical kinetic constants led to significant changes in the CO 2 production rate, while the ones in the electrochemical kinetic constants as well as in the reaction rate constants of reactions where the BBS was participating in, had minor or almost no effect on the CO 2 production rate.
Reaction rate curves have been generated to study the effect of various operating parameters on the catalytic performance. It was observed that rise in the P in CO resulted in "volcano-type" behaviour in the CO 2 NEMCA production rate, while "S-type" behaviour was favoured in the CO 2 Faradaic production rate. At low P in CO values, a sharp increase was observed in both Faradaic and Non-Faradaic rates. On the other hand, increasing the P in O 2 resulted in "S-type" behaviour in the CO 2 production rate due to both NEMCA and Faradaic effects. Moreover, the performance of the solid oxide single pellet was improved when increases in temperature were introduced in the case of variable P in O 2 for P in O 2 /P in CO o6 and the opposite effect was observed in the case of variable P in CO for higher than 1.5 kPa P in CO values. We believe that such detailed macroscopic models, aided by a good range of experiments, can greatly enhance our understanding of electrochemically promoted systems.

C i
concentration of species i in the gas phase (mol m À 3 ) D i diffusion coefficient of species i (m 2 s À 1 ) E A activation energy of a reaction (J mol À 1 ) E eq equilibrium potential difference at electrochemical reaction interface (V) F faraday constant (A s mol À 1 ) F d feed of gas mixture (m 3 s À 1 ) J i current density (A m À 2 ) J A/C anodic/cathodic current density (A m À 2 ) J 0 exchange current density (A m À 2 ) k i rate constant for adsorption (m 3 mol À 1 s À 1 ) Rate constant for surface reaction (s À 1 ) k -i rate constant for desorption (s À 1 ) k o,i pre-exponential coefficient of adsorption rate constant (m 3 mol À 1 s À 1 ) k o,-i pre-exponential coefficient of desorption rate constant (s À 1 ) L length of electrolyte/anode, m) L count length of the counter electrode in the cathode (m) M i molecular weight of species i (kg mol À 1 ) N i total molecular flux of species i (mol m À 2 s À 1 ) N s concentration of active sites on catalytic surface, mol m À 2 ) n wxp number of experimental points P i partial pressure of species i (atm) P ref reference pressure, atm) Q j charge source term (A m À 3 ) R ideal gas constant (J mol À 1 K À 1 ) R i volumetric reaction rate of species i (mol s À 1 m À 3 ) r i rate of reaction i (s À 1 ) R i overall reaction rate of species i (s À 1 ) S i sticking coefficient of species (i) T absolute temperature (K) t time (s) V OC thermodynamic open circuit potential (V) W width of electrolyte/anode/cathode (m) Greek symbols α A anode charge transfer coefficient α C cathode charge transfer coefficient γ pre-exponential coefficient (A m À 2 ) ε porosity η overpotential (V) θ i coverage of species i on Pt surface μ i chemical potential of species i (J mol À 1 ) π mathematical constant ρ charge density (A m À 3 ) s charge conductivity (Ω À 1 m À 1 ) Φ local electrostatic potential (V) Φ pellet operating potential in the solid oxide single pellet (V)