Numerical Modeling of CO₂ Reduction Reaction in a Batch-Cell with Different Working Electrodes

1. Introduction

In the recent decades, an increasing use of fossils have raised the level of greenhouse gases (GHG), especially carbon dioxide (CO₂), in the atmosphere. CO₂ concentration in atmosphere has increased from 280 ppm, from the inception of an industrial era, to an alarming 410 ppm [1,2,3]. Worldwide efforts are being devoted to curb excess amount of CO₂ from atmosphere to minimize its detrimental effect such as global warming. Carbon capture utilization and sequestration (CCUS) is a process to achieve net zero carbon emission either by CO₂ sequestration or by utilizing CO₂ as a raw feed stock to produce chemical and fuels [4]. Electrochemical conversion of CO₂ into value added chemicals is emerging as a potential technology [5]. Electrochemical reactor or electrolyzer is used to study electrochemical CO₂ reduction reaction (eCO₂RR), and some popular electrolyzers include a batch-cell, H-cell (with and without flow), flow electrolyzer (with catholyte), zero-gap membrane electrode assembly (MEA), and microfluidic reactor [6,7,8,9]. Despite substantial research in the field of electrochemistry, there are some bottlenecks such as CO₂ carbonation [10,11,12,13] and gas diffusion layer (GDL) flooding [14,15,16,17] that are preventing its realization to industrial scale. Besides a good catalyst, a robust electrolyzer design is required for optimal cell performance.

Electrochemical reduction of CO₂ is a complex physical phenomenon. It includes simultaneous physical and chemical processes such as hydrodynamics in free and porous media, mass transfer, heat transfer, gas/liquid two phase flow, gas/liquid inter-phase mass transfer, chemical and electrochemical reactions, and charge conservation. A comprehensive numerical model is required to address the intricacy associated with an electrolyzer operation. An early attempt was made by Gupta et al. [18] by presenting a mathematical model for the estimation of local CO₂ concentration and local pH in a boundary layer over a planar electrode at a certain current density. Hashiba et al. [19] presented a 1-D model for the aqueous system, assuming 100% CH₄ faradaic efficiency, with 100 $\mu$m thick hydrodynamic boundary layer near the working electrode. This model was used to demonstrate the effect of electrolyte buffer capacity on CO₂ mass transport and pH profile within the boundary layer. Corpus et al. [20] coupled the continuum model, 1-D planar electrode model with associated boundary layer, with covariance matrix adaptation evolution strategy (CMA-ES) scheme to identify the kinetic rate parameters in Tafel equation without mass transfer limitation effect. Traditional method for obtaining the kinetic rate parameters is accompanied by unquantifiable errors that leads to overprediction of current density. Bui [21] developed 1-D transient boundary layer model with copper (Cu) planar electrode to study pulsed CO₂ electrolysis. This work established that pulsing dynamically changes pH and CO₂ concentration near solid Cu electrode, which favors C₂₊ products. Qui et al. [22] developed 2-D transient models for batch-cell, planar flow cell, and gas diffusion electrode (GDE) flow cells to study the temporal and spatial variations in local CO₂ concentration and its effect on product selectivity, mainly alcohols. Other than batch-cell and H-cell, numerical models for CO₂ reduction in flow-electrolyzers have been reported in the literature, as well. Whipple et al. [23] developed a 2D numerical model for microfluidic electrolyzer to evaluate catalyst under various cell operating conditions, and to study the effect of pH on reactor efficiency. Weng et al. published three numerical models for flow electrolyzers: a 1-D GDE model [24], represents the cathodic section of flow-electrolyzer, with silver (Ag) catalyst to demonstrate the superiority of GDE over planar electrode, zero-gap membrane electrode assembly (MEA) model [25] with Ag catalyst to study gas-fed (full-MEA) and liquid-fed (exchange-MEA) configurations, and this model was extended to include copper (Cu) kinetics [26]. Kas et al. [27] developed the 2-D numerical model for GDE to emphases on the concentration gradients in the flow channels. Ehlinger et al. [28] used 1-D planar electrode model to fit the kinetics for hydrogen formation from bicarbonates; subsequently, kinetic parameters from planar electrode model were simulated in 1-D MEA model to study flooding and to perform sensitivity analysis of MEA design and operating parameters. Besides, in numerous studies numerical modeling is used to support experimental findings such as carbonation, water management, local pH and species concentration, and electrolyzer optimization [14,29,30,31,32,33,34,35].

Aforementioned references in the literature are about planar electrode (solid electrode) in batch-cell or H-cell, or flow type H-cell, and GDE in H-cell with flow configuration or in flow electrolyzers including microfluidic cell, flow electrolyzer with catholyte and zero-gap MEA. Despite the limited utility of batch-cell in the field of electrochemistry, its significance in the characterization of electrocatalyst cannot be forsaken/undermined, however, a limited number of modeling and optimization studies are available on batch-cell compared to the wealth of knowledge on flow electrolyzers. To date, there is no study available in literature on the numerical modeling of batch-cell with working electrodes other than planar electrode, whereas, most of the researched electrocatalysts are available in nano-particle form [5,36,37] that requires some support in order to be tested in an electrolyzer. Gas diffusion layer electrode (GDLE), and glassy carbon electrode (GCE) are employed to support electrocatalyst available in nano-particle form [38]. Generally, batch-cell has three electrodes including working electrode (WE) where CO₂ reduction takes place, counter electrode (CE) completes the circuit, and reference electrode (RE) that provides control over WE. Batch-cell modeling with different working electrodes could provide useful insight regarding local species concentration that influences the rate of eCO₂RR.

In this study, a 1-D, steady state, isothermal, single-phase numerical model is developed for batch-cell with three different WEs including solid electrode (SE), glassy carbon electrode (GCE), and gas diffusion layer electrode (GDLE) [39,40]. Numerical model is validated with the experimental results of CuAg_0.5Ce_0.2 catalyst, developed indigenously in our lab, in a batch-cell with GCE. Numerical model is further tested with the published experimental results of Ebaid et al. [41] for CO₂ reduction over Cu foil in a gas-tight electrolyzer with anion exchange membrane separating cathode and anode. Experimental setup details of the former catalyst including catalyst synthesis, batch-cell operation, and results analysis are reported in the supplementary information. Numerical model is then used to investigate total current density (TCD) behavior with variation in batch-cell parameters including boundary layer thickness, CL porosity, CL thickness, nature of an electrolyte and electrolyte strength. This study provides information to the experimentalist regarding the selection of working electrode, and the critical analysis of key parameters associated with the batch-cell that could be adjusted for desired performance.

Although, numerical model demonstrates the complex physical and chemical processes inside the boundary layer and the porous catalyst layer (GCE, GDLE), and it shows good agreement with experimental findings, but some assumptions have been adopted in the model that should be replaced in future work with relevant physics to capture even more realistic results. In batch-cell experiments, CO₂ saturated electrolyte is the only source of CO₂, and CO₂ concentration decreases with time, while in the model CO₂ concentration is fixed at the interface of bulk of an electrolyte and boundary layer. Boundary layer thickness is a strong function of stirring speed, higher stirring speed results in a thinner boundary layer and vice versa is possible, but in the model boundary layer thickness is adopted from the literature [7,21,24,27,42], and the impact of various boundary layer thickness have been discussed in detail. Cu based electrocatalysts are known for producing range of gaseous and liquid CO₂ reduction products. During the batch-cell operation, there is a high probability that some of the CO₂ reduction products interact with each other, and influence pH, CO₂ solubility, and species transport. To avoid the model complexity, it is assumed that liquid products do not participate in any reaction. Temperature variations are common in a batch-cell mainly due to the application of applied voltage and heat of reactions, and it could impact the rate of eCO₂RR, however, isothermal conditions have been assumed for this study. This discussion highlights potential pathways in batch-cell modeling research work.

2. Numerical Model

The region of interest for numerical modeling in batch-cell is a WE (glassy carbon electrode), where CO₂ is reduced. GCE has a solid flat surface, coated with CL. Boundary layer (BL) is formed on the surface of catalyst that governs mass transfer of species between bulk of an electrolyte and WE, while bulk of an electrolyte achieves uniformity due to constant stirring. Therefore, computational domain (CD) for GCE has two important regions: porous catalyst layer (CL), and BL. A schematic of CD is shown in the Figure 1 with all the necessary boundary conditions. The CD is modeled as 1D.

A comprehensive numerical model is developed to deconvolute the physio-electrochemical phenomenon such as mass transfer due to diffusion and electromigration, charge conservation, electrochemical reaction, and chemical reactions in the CD. A total of 6 species have been modeled: dissolved CO₂, hydroxyl ions (OH^-1), bicarbonate ions (HCO₃^1-), carbonate ions (CO₃^2-), protons (H⁺), and potassium ions (K⁺). Catalyst used in this study reduces CO₂ to five different products with higher selectivity for C₂₊ products: carbon mono-oxide (CO), methane (CH₄), ethylene (C₂H₄), ethanol (C₂H₅OH), and propanol (C₃H₇OH). Alongside these products, a competing parasitic hydrogen evolution reaction (HER) also takes place. The electrochemical reaction for each of these products is given in the Table 1 [36]. It is assumed that gaseous products (H₂, CO, CH₄, C₂H₄) bubble out as soon as they are produced, while liquid phase products act as neutral species, and do not participate in any reaction. Water (H₂O) is assumed to be ubiquitous specie in CD; therefore, it is not included in the modeling. Electrolyte is a 1 M KOH solution. Besides electrochemical conversion of CO₂ into useful products, a considerable amount of CO₂ reacts reversibly with the hydroxyl ions (OH^-1) to produce bicarbonates (HCO₃^-1), and carbonates (CO₃^2-) ions. These reactions govern local pH in CD. Due to high pH (~1 M KOH pH=14) of an electrolyte, only the reactions in the Table 2 have been considered in the model [43]. Precipitation of bicarbonate (KHCO₃), and carbonate (K₂CO₃) salts are ignored in the numerical model. Although, salt build-up has a substantial deteriorating effect on the CL [11,12] and it also increases required cell potential if experiment is carried out for extended period [10], but in case of batch-cell modeling, this phenomenon can be ignored as a fresh electrolyte solution is used for each batch cell experiment [12].

Specie balance is used to model concentration of species in the CL, and the BL.

$$\nabla N_j=R_{c,j}+R_{e,j}$$

(1)

N_j describes the molar flux of specie j in CD, while source terms on the right hand represents the consumption and the generation of species due to chemical and electrochemical reactions, respectively. Nernst-Planck equation is invoked for the molar flux of species in CD [40].

$$N_j=-D_j\nabla c_j+{\displaystyle \frac{z_{j}F}{RT}}{D}_j{c}_j\nabla {\varphi }_L$$

(2)

First term on the right-hand side is the diffusion of species according to the Fick’s law, while the second term represents the electro-migration of charged species; therefore, the second becomes zero for neutral species such as dissolved CO₂. The diffusion of species in porous CL is corrected using Bruggeman correlation [44]. The effective diffusion coefficient ($D_j^{eff}$) of species in porous CL depends on its porosity (${\epsilon }_{CL}$), and tortuosity (${\tau }_{CL}$).

$$D_j^{eff}=D_j{\displaystyle \frac{{\epsilon }_{CL}}{{\tau }_{CL}}}$$

(3)

$${\tau }_{CL}={\displaystyle \frac{1}{\sqrt{{\epsilon }_{CL}}}}$$

(4)

An additional equation is required for electrolyte potential (${\varphi }_l$) to close the degree of freedom; therefore, electroneutrality equation is included in the model.

$$\sum _j{z}_j{c}_j=0$$

(5)

Nernst-Planck equation assumes dilute-solution theory [24,40]. Modeled domain is not dilute under specified conditions; however, it is a reasonably good assumption as the TCD is around 100 mA cm^-2. Concentrated solution theory requires additional diffusion coefficients which are not available in the literature.

The chemical consumption/generation term is defined as follows;

$$R_{c,j}=\sum _n{v}_{j,n}\left(k_{f,n}\prod _{v_{j,n<0}}{c}_j^{-v_{j,n}}-k_{r,n}\prod _{v_{j,n>0}}{c}_j^{v_{j,n}}\right)$$

(6)

$$K_n={\displaystyle \frac{k_{f,n}}{k_{r,n}}}$$

(7)

$v_{j,n}$ is the stoichiometric coefficient. It is negative for reactant species, and positive for product species. $k_{f,n}$ and $k_{r,n}$ are the forward and the reverse reaction rate constants, respectively. $K_n$ is an equilibrium constant.

The electrochemical source term in specie balance is defined as;

$$R_{e,j}=-M_j\sum _k{\displaystyle \frac{v_{j,k}{a}_s{i}_k}{n_{k}F}}$$

(8)

The term ($M_j$) represents molar weight of species. $a_s$ is the effective surface area available for CO₂ reduction in porous CL domain, $n_k$ represents the electron transferred in each reaction $k$, $i_k$ is the partial current density, and F is the Faraday’s constant. The electrode kinetics or electrochemical reaction kinetics are modeled using concentration dependent Butler-Volmer (BV) equation [45]. The BV equation describes both anodic half reaction, and cathodic half reaction. The BV equation reduces to Tafel equation that describes the electrochemical conversion of dissolved CO₂ into products at WE. A detailed derivation of the Tafel equation is provided in the supplementary information. The final form of the Tafel kinetic equation used in the model is;

$$i_k=-i_{0,k}{\left({\displaystyle \frac{c_{C{O}_2}}{c_{ref}}}\right)}^{{\gamma }_{C{O}_2,k}}\exp \left(-{\gamma }_{pH,k,SHE}pH\right)\exp \left[-{\displaystyle \frac{{\alpha }_{c,k}F}{RT}}\left({\varphi }_s-{\varphi }_l-E_{0,k}\right)\right]$$

(9)

This formulation of Tafel kinetic equation is adopted from Bui [21]. Tafel equation has several unknown parameters that are either the direct inputs to simulation or obtained through data fit to experimental findings. ${\varphi }_s$ is applied cell voltage, while ${\varphi }_l$ is evaluated by solving Nernst-Planck equation, and electroneutrality equation (see equation 2, and 5). $E_{0,k}$ is the thermodynamic potential available in literature for almost all product species. ${\gamma }_{C{O}_2}$ represents the reaction order of CO₂ for specific product specie. It could be obtained either through experimental data fitting of CO₂ reduction at various applied potential as reported in previous studies by Weng [25] or through rigorous microkinetic modeling discussed by Rae [42]. According to the density functional theory (DFT) calculations, ${\gamma }_{C{O}_2}$ is the number of elementary steps in reduction of CO₂ towards certain product specie [21], which is the most simplified way to obtain this parameter. In this model, ${\gamma }_{C{O}_2}$ for all product species are taken from the reference [21], and it is listed for each CO₂ reduction product in Table 4. ${\gamma }_{pH}$ is the pH rate ordering parameter. It is evaluated through experimental data fit of current density vs. pH at a specified applied voltage [21]. ${\gamma }_{pH}$ values for each CO₂ reduction product are listed in Table 4. $c_{C{O}_2}$ and $pH$ are local CO₂ concentration, and local pH, respectively. These values are quite challenging to obtain experimentally; therefore, simulation is performed using experimental input to get these values. $c_{ref}$ is the reference concentration that is set to 1 M. ${\alpha }_{c,k}$ is the transfer coefficient, while $i_{0,k}$ is the exchange current density. More on these parameters will be discussed in the next section.

Charge conservation in CD is modeled using the following equations.

$$\nabla ·i=0$$

(10)

$$\nabla ·i_s=-\nabla ·i_l=S$$

(11)

$$S=-a_s\sum _k{i}_k$$

(12)

Ohm’s law is used to define the relation between current and voltage.

$$i={\sigma }_m\nabla \mathsf{\varphi }$$

(13)

${\sigma }_m$ reparents the conductivity solid catalyst, and liquid electrolyte. In porous CL, effective conductivity is used according to the Bruggeman correlation.

$${\sigma }^{eff}=\sigma {\displaystyle \frac{{\epsilon }_{CL}}{{\tau }_{CL}}}$$

(14)

$${\tau }_{CL}={\displaystyle \frac{1}{\sqrt{{\epsilon }_{CL}}}}$$

(15)

Concentration of dissolved CO₂ in 1 M KOH at STP is 24.57 mM. It is estimated using Henry’s law, and it is corrected for electrolyte strength using Sechenov’s constants [46].

$$C_{0,c{o}_2}^{aq}=H_{C{O}_2}{C}_{0,c{o}_2}^{gas}$$

(16)

$H_{C{O}_2}$ is the Henry’s constant for CO₂, which is defined as [47,48];

$$\ln \left(H_{C{O}_2}\right)=93.4517\times {\displaystyle \frac{100}{T}}-60.2409+23.3585\times \ln \left({\displaystyle \frac{T}{100}}\right)$$

(17)

In case of an electrolyte with high concentration of ions, solubility of CO₂ is reduced. Reduced CO₂ solubility in an electrolytic solution is evaluated using Sechenov’s constants [46].

$$\log \left({\displaystyle \frac{C_{0,c{o}_2}^{aq}}{C_{c{o}_2}^{aq}}}\right)=K_s{C}_s$$

(18)

$K_s$ is the Sechenov’s constant, and $C_s$ is the molar concentration of ions. Sechenov’s constants value are provided in Table 3.

$$K_s=\sum \left(h_{ion}+h_G\right)$$

(19)

$$h_G=h_{G,0}+h_T\left(T-298.15\right)$$

(20)

At x=L_CD, the concentration of species is set to the bulk of an electrolyte. CO₂ starts diffusing from the interface of BL and bulk of an electrolyte to reach the surface of CL, where it undergoes eCO₂RR. Since, CL is a porous domain, it offers more active sites for CO₂ reduction, but due to intricate network of pores, diffusion of CO₂ inside CL is reduced. The nature of eCO₂RR products depend on the catalyst. The thickness/length of BL is about 100 $\mu m$ [21]. However, this value changes with the stirring speed (in RPM). Voltage (${\varphi }_s=E_{WE}$) is applied at x=0, whereas it is assumed that RE is located at the interface of BL and bulk of an electrolyte (x=L_CD), where electrolyte potential is set to zero (${\varphi }_l=0$). Moreover, steady state, and isothermal conditions have been assumed.

2.1. Kinetic Parameters from Experimental Data

In this section, Tafel equation parameters for each CO₂ reduction product are obtained from the experimental data. Experimental data is provided in Table S1 (refer to supplementary information). Tafel kinetic equation discussed in the last section (see equation 9) is used to estimate kinetic parameters. Equation 9 has six unknown parameters including ${\gamma }_{C{O}_2,k}$, ${\gamma }_{pH,k,SHE}$, $c_{C{O}_2}$, $pH$, $i_{0,k}$, ${\alpha }_{c,k}$. Parameter values for ${\gamma }_{C{O}_2,k}$, and ${\gamma }_{pH,k,SHE}$ are provided in Table 4. Parameters $c_{C{O}_2}$, and $pH$ are listed in Table 5. $i_{0,k}$, and ${\alpha }_{c,k}$ values are reported in Table 6.

Local CO₂ concentration, and pH values are estimated by simulating experimental data. Experimental TCD is set as a boundary condition. Mass flux boundary condition, estimated from experimental partial current density (PCD), is used for dissolved CO₂ and hydroxyl ions (OH^-1) species balance. In this way, local average CO₂ concentration, and local average pH values are estimated in CD. These values are provided in Table 5.

Remaining two parameters $i_{0,k}$, and ${\alpha }_{c,k}$ are obtained through comparing with experimental data. Information provided in Table 4, and Table 5 is used to evaluate the normalized current density.

$$i_k^{norm}=-{\displaystyle \frac{i_k^{exp}}{\exp \left(-{\gamma }_{pH,k,SHE}pH\right){\left(\frac{c_{C{O}_2}}{c_{ref}}\right)}^{{\gamma }_{C{O}_2,k}}}}$$

(21)

Tafel equation is simplified as;

$$i_k^{norm}=i_{0,k}\exp \left(-{\displaystyle \frac{{\alpha }_{c,k}F}{RT}}\eta \right)$$

(22)

$$\ln \left(i_k^{norm}\right)=\ln \left(i_{0,k}\right)+\left(-{\displaystyle \frac{{\alpha }_{c,k}F}{RT}}\right)\eta$$

(23)

The slope, and the intercept of the equation 23 provide ${\alpha }_{c,k}$, and $i_{0,k}$, respectively. Normalized current density plots as a function of overpotential for each experimentally recorded CO₂ reduction product are shown in Figure S3.

2.2. Numerical Model for Solid Electrode (SE) and Gas Diffusion Layer Electrode (GDLE)

In this section, numerical model for other two WEs have been presented: solid electrode (SE), and gas diffusion layer electrode (GDLE). In practical, SE is made up of a single metal such as copper (Cu), but, in this study, it is assumed that SE is made-up of catalyst CuAg_0.5Ce_0.2 used in the experiments. GDLE is a GDL paper that contains both macro pores (carbon fibers), and micro pores (carbon powder). A layer of catalyst is coated on top of micro pores layer of GDL paper. Electrochemical kinetic parameters obtained in the previous section are used to simulate both electrodes.

2.2.1. Solid Electrode Numerical Model

A schematic of CD for SE is shown in Figure 2. There is only boundary layer (BL) in the CD of SE. CD contains six species: dissolved CO₂, OH^-, HCO₃^-, CO₃⁼, H⁺, K⁺. Concentration of these species is set to the bulk of an electrolyte at (x=L_CD). CO₂ diffuses from the interface of BL and the bulk of an electrolyte (x=L_CD) to reach the surface of the electrode. Electrochemical reaction happens at the surface of the electrode (x=0). Therefore, Neumann flux boundary condition is used to model the consumption and production of CO₂ and OH^- ions at the surface of an electrode (x=0), respectively. Voltage (${\varphi }_s$) is applied at surface of the electrode, while electrolyte potential (${\varphi }_l$) is set to 0 at x=L_CD, assuming that RE is placed at this location. In case of SE, length of the CD is similar to the thickness of the BL, and the length of the CD is 100 $\mu$m. All the other parameters are similar to the GCE.

2.2.2. Gas Diffusion Layer Electrode Numerical Model

GDLE has a porous structure. Gas diffusion layer (GDL) paper has a layer of macro pores (carbon fibers), and a thin layer of micro pores (carbon power). However, it is assumed that GDL paper has a uniform structure with same pore size, and it has isotropic porosity. CL that is coated on top of GDLE has a porous structure, as well. GDLE computational domain is shown in Figure 3. Porous nature of GDLE allows reactant species to diffuse from both front (with catalyst layer), and back (without catalyst layer). CD has 4 regions: two boundary layers, porous GDL, and porous CL. CD contains dissolved CO₂, OH^-, HCO₃^-, CO₃⁼, H⁺, and K⁺. Concentration of species is set to bulk on an electrolyte at x=0, and x=L_CD. CO₂ diffuses from both ends of the CD to reach CL, where it undergoes eCO₂RR. GDL serves two important functions: it supports the CL, and acts as a medium for electronic conductor. Voltage (${\varphi }_s$) is applied at the interface of GDL/BL. It is assumed that RE is located at x=L_CD, and the electrolyte potential (${\varphi }_l$) at this location is set to zero.

2.3. Numerical Simulation

Numerical model developed in this study is simulated using COMSOL Multiphysics v6.1 environment. COMSOL uses finite element technique to discretize the differential equations. MUMPS direct general solver is used to solve set of algebraic equations (Ax=B). Total number of nodes in CD are 200; number of nodes in CL domain, and BL domain are set to 150, and 50, respectively. Information required to simulate the model is provided in Table 6.

Grid independent study is performed to check the model consistency and robustness. CL domain is more sensitive to species concentration gradient than the BL domain, as eCO₂RR happens inside the pores of CL; therefore, smaller element size (0.30 $\mu$m) is used CL domain, comparing with BL domain with relatively bigger element size (0.67 $\mu$m). Grid independent solution is obtained when the variation in the TCD, at a fixed applied cathodic voltage of -1.0 V vs RHE, goes below 1% by increasing number of nodes (decreasing element size) in CD. Grid independent solution is shown in Figure 4.

2.4. Numerical Model Validation

Numerical model is validated by comparing simulation TCD result with experimental data, see Figure 5. Numerical model shows a good agreement with experiment, however there is a slight difference in simulated TCD slope and experimental TCD slope. Numerical model predicts TCD with positive slope, however, the slope of experimental data shows a decreasing trend with the increase in applied cathodic voltage. In experiments, electrolyte saturated with CO₂ is used, and batch-cell is sealed for the duration of an experiment. Concentration of dissolved CO₂ in an electrolyte becomes a function of time, and the concentration of CO₂, which is 24.57 mM in 1 M KOH at STP, keeps on decreasing. This physical phenomenon is not catered in the simulation, and it is simply replaced with a Dirichlet concentration boundary condition. This assumption maintains a relatively higher, and constant concentration gradient for CO₂ at a certain applied cathodic voltage, than the experiments where mass transfer resistance increases with decreasing CO₂ concentration in the bulk of an electrolyte. Therefore, decreasing slope trend for experimental current density is due to scarce supply of CO₂, which becomes significant at higher applied voltages. Furthermore, use of raw technique for kinetic parameter estimation, and incapability of 1-D model to fully describe underlying physics could also affect the numerical simulation results. Despite the simplification of the numerical model, it is able to predict TCD with quite a good accuracy.

A comparison of simulated PCD of eCO₂RR products with experimental findings is shown in Figure S4. Except for ethylene, all the individual current density shows a good match with experiment. Experimental data for ethylene shows irregular behavior, especially the data point at -0.7 V vs RHE and -1.0 V vs RHE, which are significantly away from the normal trend.

Numerical model developed in this study is also validated with the experimental findings of Ebaid et al [41]. CU foil and graphite rod was used as a working electrode and counter electrode, respectively. Electrolyte was 0.1 M CsHCO₃ saturated with CO₂. Validation results for TCD and PCD are shown in Figure 6 and Figure S5, respectively.

Table 6. Batch-cell model parameters.

Parameter		Value	Unit		Reference
Design Parameters
LCL		15	$\mu m$		Experimental input
LBL		100	$\mu m$		Fitted
LGDL		325	$\mu m$		Experimental input
T		298	K		This study
EWE		-0.7 – -1.1	V vs RHE
Electrochemical Kinetic Rate Parameters for CuAg0.5Ce0.2 Catalyst
$i_{0,h2}$		97.3	mA cm-2		This Study
$i_{0,CO}$		8.89 x 108	mA cm-2
$i_{0,C{H}_4}$		9.50 x 105	mA cm-2
$i_{0,C_2{H}_4}$		2.8	mA cm-2
$i_{0,EtOH}$		4.10 x 10-2	mA cm-2
$i_{0,PrOH}$		9.50 x 10-3	mA cm-2
${\alpha }_{H_2}$		0.0274
${\alpha }_{CO}$		0.1001
${\alpha }_{C{H}_4}$		0.1391
${\alpha }_{C_2{H}_4}$		0.1222
${\alpha }_{EtOH}$		0.1490
${\alpha }_{PrOH}$		0.1676
$E_{0,H_2}$		0	V vs RHE		[36]
$E_{0,CO}$		-0.11	V vs RHE
$E_{0,{CH}_4}$		0.17	V vs RHE
$E_{0,{C_{2}H}_4}$		0.07	V vs RHE
$E_{0,EtOH}$		0.08	V vs RHE
$E_{0,PrOH}$		0.09	V vs RHE
Chemical Reaction Rate Parameters
$k_{1,f}$		5.93 x 10-3	m3 mol-1s-1		[27]
$k_{2,f}$		1.0 x 10-8	m3 mol-1s-1		[27]
$k_{3,f}$		0.04	s-1		[24]
$k_{4,f}$		56.281	s-1		[24]
$k_{w,f}$		0.001	mol L-1s-1		[26]
$k_{1,b}$		1.34 x 10-4	s-1		[27]
$k_{2,b}$		2.15 x 10-4	s-1		[27]
$k_{3,b}$		9.37 x 104	L mol-1s-1		[24]
$k_{4,b}$		1.23 x 1012	L mol-1s-1		[24]
$k_{w,b}$		1.0 x 1011	L mol-1s-1		[26]
Porous Media Properties
${\epsilon }_{GDL}$		0.72			This Study
${\epsilon }_{CL}$		0.3			This Study
${\sigma }_{GDL}$		220	S m-1		[24]
${\sigma }_{CL}$		100	S m-1		[24]
$a_s$		1.0 x 107	m-1		This Study
Liquid Species Diffusion Coefficients
$D_{H^{+}}$	$4.49\times {10}^{-9}\exp \left(-1430\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{273.15}}\right)\right)$			m2 s-1		[49]
$D_{O{H}^{-}}$	$2.89\times {10}^{-9}\exp \left(-1750\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{273.15}}\right)\right)$			m2 s-1		[49]
$D_{HC{O}_3^{-}}$	$7.016\times {10}^{-9}{\left({\displaystyle \frac{T}{204.0282}}-1\right)}^{2.3942}$			m2 s-1		[50]
$D_{C{O}_3^{=}}$	$5.447\times {10}^{-9}{\left({\displaystyle \frac{T}{210.2646}}-1\right)}^{2.1929}$			m2 s-1		[50]
$D_{C{O}_2}$	$2.17\times {10}^{-9}\exp \left(-2345\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{303}}\right)\right)$			m2 s-1		[51]
$D_{K^{+}}$	$1.957\times {10}^{-9}\exp \left(-2300\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{298.15}}\right)\right)$			m2 s-1		[51]

Table 6. Batch-cell model parameters.

Parameter		Value	Unit		Reference
Design Parameters
LCL		15	$\mu m$		Experimental input
LBL		100	$\mu m$		Fitted
LGDL		325	$\mu m$		Experimental input
T		298	K		This study
EWE		-0.7 – -1.1	V vs RHE
Electrochemical Kinetic Rate Parameters for CuAg0.5Ce0.2 Catalyst
$i_{0,h2}$		97.3	mA cm-2		This Study
$i_{0,CO}$		8.89 x 108	mA cm-2
$i_{0,C{H}_4}$		9.50 x 105	mA cm-2
$i_{0,C_2{H}_4}$		2.8	mA cm-2
$i_{0,EtOH}$		4.10 x 10-2	mA cm-2
$i_{0,PrOH}$		9.50 x 10-3	mA cm-2
${\alpha }_{H_2}$		0.0274
${\alpha }_{CO}$		0.1001
${\alpha }_{C{H}_4}$		0.1391
${\alpha }_{C_2{H}_4}$		0.1222
${\alpha }_{EtOH}$		0.1490
${\alpha }_{PrOH}$		0.1676
$E_{0,H_2}$		0	V vs RHE		[36]
$E_{0,CO}$		-0.11	V vs RHE
$E_{0,{CH}_4}$		0.17	V vs RHE
$E_{0,{C_{2}H}_4}$		0.07	V vs RHE
$E_{0,EtOH}$		0.08	V vs RHE
$E_{0,PrOH}$		0.09	V vs RHE
Chemical Reaction Rate Parameters
$k_{1,f}$		5.93 x 10-3	m3 mol-1s-1		[27]
$k_{2,f}$		1.0 x 10-8	m3 mol-1s-1		[27]
$k_{3,f}$		0.04	s-1		[24]
$k_{4,f}$		56.281	s-1		[24]
$k_{w,f}$		0.001	mol L-1s-1		[26]
$k_{1,b}$		1.34 x 10-4	s-1		[27]
$k_{2,b}$		2.15 x 10-4	s-1		[27]
$k_{3,b}$		9.37 x 104	L mol-1s-1		[24]
$k_{4,b}$		1.23 x 1012	L mol-1s-1		[24]
$k_{w,b}$		1.0 x 1011	L mol-1s-1		[26]
Porous Media Properties
${\epsilon }_{GDL}$		0.72			This Study
${\epsilon }_{CL}$		0.3			This Study
${\sigma }_{GDL}$		220	S m-1		[24]
${\sigma }_{CL}$		100	S m-1		[24]
$a_s$		1.0 x 107	m-1		This Study
Liquid Species Diffusion Coefficients
$D_{H^{+}}$	$4.49\times {10}^{-9}\exp \left(-1430\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{273.15}}\right)\right)$			m2 s-1		[49]
$D_{O{H}^{-}}$	$2.89\times {10}^{-9}\exp \left(-1750\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{273.15}}\right)\right)$			m2 s-1		[49]
$D_{HC{O}_3^{-}}$	$7.016\times {10}^{-9}{\left({\displaystyle \frac{T}{204.0282}}-1\right)}^{2.3942}$			m2 s-1		[50]
$D_{C{O}_3^{=}}$	$5.447\times {10}^{-9}{\left({\displaystyle \frac{T}{210.2646}}-1\right)}^{2.1929}$			m2 s-1		[50]
$D_{C{O}_2}$	$2.17\times {10}^{-9}\exp \left(-2345\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{303}}\right)\right)$			m2 s-1		[51]
$D_{K^{+}}$	$1.957\times {10}^{-9}\exp \left(-2300\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{298.15}}\right)\right)$			m2 s-1		[51]

4. Conclusions

In conclusion, we have developed a numerical model to study CO₂ reduction reaction in a three electrodes batch-cell with three different working electrodes (glassy carbon electrode, solid electrode, gas diffusion layer electrode). Numerical model is simulated using COMSOL software. Following are key conclusion from our numerical simulation.

• Thickness of boundary layer (BL) has an inverse relation with the total current density (TCD), and hence with the performance of batch-cell. Thicker BL offers more resistance to CO₂ diffusion, while thinner BL improves the mass transfer of other species such as hydroxyl ions (OH^-). Magnetic stirrer speed could be manipulated to control the BL thickness in a batch-cell.
• Catalyst layer (CL) porosity also demonstrates an inverse relation with TCD; higher the CL porosity, lower the TCD will be, and vice versa. However, higher CL porosity reduces the mass transfer resistance. An optimum selection of CL porosity is required to balance the effective diffusion rate with good TCD.
• CL thickness has a significant effect on the performance of the batch-cell. A thinner CL is much more effective in terms of overall batch cell performance, compared with a relatively thicker CL. Batch-cell suffers from poor CO₂ diffusion rate, and limited amount of CO₂ is available for electrochemical CO₂ reduction reaction (eCO₂RR), therefore, only a small portion of CL is used in eCO₂RR. A thinner CL saves catalyst, improves species mass transfer rate, and results in a higher current density. An optimum balance of catalyst loading, CL thickness, and CL porosity could enhance the overall cell performance.
• Electrolyte nature i-e cations and anions impacts the CO₂ reduction reaction. 1 M KOH performs well in comparison with 1 M KHCO₃ in a batch-cell. KHCO₃ is a buffer solution, therefore, it consumes relatively more CO₂ to maintain its buffer. Batch-cell becomes deficient in CO₂, and as a result, the rate of parasitic hydrogen evolution reaction (HER) increases substantially. KHCO₃ has a near neutral pH, therefore, it will promote HER, eventually. Hence, KOH is an optimum choice, as it provides a strong alkaline environment that mitigates the probability of HER.
• KOH electrolyte with higher molarity increases the TCD define, by minimizing the charge transfer resistance between the electrodes. Nevertheless, high molarity of KOH increases the chemical conversion of CO₂ into bi-carbonates, and carbonates. An optimum distance between the electrodes could eliminate the requirement of high molar electrolyte for conductivity improvement.
• Among the three electrodes, GDLE and GCE out-performs SE in terms of TCD. In particular, the porous flow-through structure of GDLE facilitates more CO₂ into the CL, and the rate of eCO₂RR increases.