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

1. Introduction

In recent decades, an increasing use of fossils has raised the level of greenhouse gases, especially carbon dioxide (CO₂), in the atmosphere. The atmospheric CO₂ concentration has risen dramatically from 280 ppm to a concerning 410 ppm since the onset of the industrial era [1,2,3]. Therefore, a pathway to a more sustainable society is essential [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 (WE). 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 1-D continuum planar electrode model with covariance matrix adaptation evolution strategy 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 et al. [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 2-D 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 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 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] employed a 1-D planar electrode model to determine the kinetics of hydrogen formation from bicarbonates. The kinetic parameters obtained from the planar electrode model were then implemented in a 1-D MEA model to investigate flooding and conduct a 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,36,37].

Numerical modeling of batch-cells with various WEs can offer valuable insights into local species concentrations that affect the rate of eCO₂RR. To date, there is no study available in literature on the numerical modeling of batch-cell with WEs other than planar electrode, whereas most of the researched electrocatalysts are available in nano-particle form [5,38,39] 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 [40]. Aforementioned references in the literature are about planar electrode (solid electrode) in batch-cell, H-cell, and flow type H-cell, and GDE in H-cell with flow configuration, in microfluidic cell, in flow electrolyzer with catholyte and in zero-gap MEA. Although batch-cells have limited applications in electrochemistry, their importance in electrocatalyst characterization remains indispensable. However, compared to the extensive research available on flow electrolyzers, studies focusing on the modeling and optimization of batch-cells are relatively scarce.

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), GCE, and GDLE [41,42]. 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. [43] 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, catalyst layer (CL) porosity, CL thickness, nature of an electrolyte and electrolyte strength. This study provides information to the experimentalist regarding the selection of WE, 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 CL (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,44], 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 potential 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 the WE of GCE, where CO₂ is reduced. GCE has a solid flat surface, coated with CL. Boundary layer 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 for GCE has two important regions: porous CL, and boundary layer. A schematic of one-dimensional (1-D) computational domain is shown in the Figure 1 with all the necessary boundary conditions.

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 computational domain. 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 [38]. 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 species in computational domain; 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 computational domain. 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 [45]. Precipitation of bicarbonate/carbonate 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 boundary layer.

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

(1)

N_j describes the molar flux of specie j in computational domain, 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 computational domain [42].

$$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 Fick’s law, while the second term represents the electro-migration of charged species; therefore, the second term becomes zero for neutral species such as dissolved CO₂. The diffusion of species in porous CL is corrected using Bruggeman correlation [46]. 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=1-6}{z}_j{c}_j=0$$

(5)

Nernst-Planck equation assumes dilute-solution theory [24,42]. Modeled domain is not dilute under specified conditions. 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=1-3}{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=1-6}{\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 (PCD), and F is the Faraday’s constant. The electrode kinetics are modeled using concentration dependent Tafel equation [47] that describes the electrochemical conversion of dissolved CO₂. 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 et al. [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 the applied cell potential, 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. ${\gamma }_{C{O}_2}$ represents the reaction order of CO₂ for specific product species. It could be obtained either through experimental data fitting of CO₂ reduction at various applied potential [25] or through rigorous microkinetic modeling [44]. According to the density functional theory calculations, ${\gamma }_{C{O}_2}$ is the number of elementary steps in CO₂ reduction towards certain product specie, which is the most simplified way to obtain this parameter. ${\gamma }_{pH}$ is the pH rate ordering parameter. It is estimated by fitting experimental data of current density versus pH at a given applied potential. In this study, ${\gamma }_{C{O}_2}$ and ${\gamma }_{pH}$ for all product species are sourced from the reference [21], and their values are listed in Table 4. $c_{C{O}_2}$ and $pH$ are local CO₂ concentration, and local pH, respectively. Obtaining these values experimentally is quite challenging; hence, simulations are conducted, similar to Gupta et al. [18], using experimental inputs to determine them. $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.

Charge conservation in computational domain is modeled using the following equations.

$$\nabla \bullet i=0$$

(10)

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

(11)

$$S=-a_s\sum _{k=1-6}{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$ represents the conductivity of 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)

Concentration of dissolved CO₂ in 1 M KOH at STP is 24.57 mM. It is estimated using Henry’s law.

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

(15)

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

$$\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)$$

(16)

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 [50].

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

(17)

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

$$K_s=\sum _{ion=1-4}\left(h_{ion}+h_G\right)$$

(18)

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

(19)

At x=L_CD, the concentration of species is set to the bulk of an electrolyte. CO₂ starts diffusing from the interface of boundary layer and bulk 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 depends on the catalyst. The thickness/length of boundary layer is about 100 $\mu m$ [21]. However, this value changes with the stirring speed (in RPM). Potential (${\varphi }_s=E_{WE}$) is applied at x=0, whereas it is assumed that reference electrode (RE) is located at the interface of boundary layer and bulk 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.

Table 4. List of CO₂ reaction rate order parameter (${\gamma }_{C{O}_2}$) and pH rate order parameter (${\gamma }_{pH}$) values. Parameter values are sourced from Bui et al. [21].

Product Species	${\mathit{\gamma }}_{\mathit{C}{\mathit{O}}_2}$	${\mathit{\gamma }}_{\mathit{p}\mathit{H}}$
H2	0	0.5
CO	1.5	1.56
CH4	0.84	1.56
C2H4	1.36	0
EtOH	0.96	0
PrOH	0.96	0

Table 4. List of CO₂ reaction rate order parameter (${\gamma }_{C{O}_2}$) and pH rate order parameter (${\gamma }_{pH}$) values. Parameter values are sourced from Bui et al. [21].

Product Species	${\mathit{\gamma }}_{\mathit{C}{\mathit{O}}_2}$	${\mathit{\gamma }}_{\mathit{p}\mathit{H}}$
H2	0	0.5
CO	1.5	1.56
CH4	0.84	1.56
C2H4	1.36	0
EtOH	0.96	0
PrOH	0.96	0

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 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 computational domain. 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}}}}$$

(20)

Tafel equation is simplified as:

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

(21)

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

(22)

The slope, and the intercept of the equation 22 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 SE and GDLE

In this section, numerical model for other two WEs have been presented: SE, and GDLE. In practical, SE is made up of a single metal such as Cu, but, in this study, it is assumed that SE is made-up of catalyst CuAg_0.5Ce_0.2. 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 pore layer. Electrochemical kinetic parameters obtained in the previous section are used to simulate both electrodes.

2.2.1. Solid Electrode Numerical Model

Figure 2 illustrates the 1-D computational domain of SE, which comprises solely the boundary layer. Computational domain contains six species: dissolved CO₂, OH^-, HCO₃^-, CO₃⁼, H⁺, K⁺. Species flux is solved using equation 1. Concentration of these species is set to the bulk electrolyte values at (x=L_CD). CO₂ diffuses from the interface of boundary layer and bulk 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. Potential (${\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 computational domain is similar to the thickness of the boundary layer, and the length of the computational domain is 100 $\mu$m.

2.2.2. GDLE Numerical Model

GDLE consists of porous GDL and porous CL. Both GDL and CL are assumed to have a uniform structure with consistent pore size and isotropic porosity. A 1-D computational domain of GDLE is shown in Figure 3. The porous nature of the GDLE enables reactant species to diffuse from both the front side (with the CL) and the back side (without the CL). GDLE computational domain has 4 regions: two boundary layers, a porous GDL, and a porous CL. Computational domain contains dissolved CO₂, OH^-, HCO₃^-, CO₃⁼, H⁺, and K⁺, and species transport is modeled using equation 1. Concentration of species is set to bulk on an electrolyte at x=0, and x=L_CD. CO₂ diffuses from both ends of the computational domain 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. Potential (${\varphi }_s$) is applied at the interface of GDL and boundary layer. 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 Multiphysics 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 computational domain are 200; number of nodes in CL domain, and boundary layer 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 boundary layer domain, as eCO₂RR happens inside the pores of CL; therefore, smaller element size (0.30 $\mu$m) is used in CL domain, comparing with boundary layer 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 potential of -1.0 V vs RHE, goes below 1% by increasing number of nodes (decreasing element size) in computational domain. Grid independent solution is shown in Figure 4.

2.4. Numerical Model Validation

Numerical model is validated by comparing simulated TCD result with experimental data, see Figure 5. Numerical model shows 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 potential. 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 potential, than the experiments where mass transfer resistance increases with decreasing CO₂ concentration in the bulk electrolyte. Therefore, decreasing slope trend for experimental current density is due to scarce supply of CO₂, which becomes significant at higher applied potentials. 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 data 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 data of Ebaid et al. [43]. Cu foil and graphite rod were used as a WE and counter electrode (CE), 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.

3. Results and Discussion

In this section, numerical model is used to study the effect of critical parameters of batch-cell with GCE using the kinetic parameters of CuAg_0.5Ce_0.2 catalyst. The thickness of the boundary layer is a crucial factor that governs the diffusion rate of species. Thicker boundary layer offers significant resistance to CO₂ diffusion. It does not only increase diffusion time of CO₂ to reach catalyst surface; the probability of CO₂ consumption by side reactions also increases, even further reducing effective concentration of CO₂ available for eCO₂RR. Additionally, at lower cathodic potentials, an adequate amount of CO₂ is available for eCO₂RR. However, the TCD sharply decreases at higher cathodic potentials due to an imbalance between the rate of eCO₂RR and CO₂ mass transfer. The effect of boundary layer thickness on TCD at -1.0 V vs RHE is shown in the Figure 7. TCD is obtained for 4 different boundary layer thicknesses of equal interval: 40 $\mu$m, 80 $\mu$m, 120 $\mu$m, 160 $\mu$m. Figure 7 shows an inverse relationship between TCD, and boundary layer thickness, and the trend is not linear. A substantial improvement in the TCD at boundary layer thickness of 40 $\mu$m, comparing with 160 $\mu$m, is due to the lower resistance to CO₂ mass transfer and decline in chemical conversion of CO₂ into bi-carbonates and carbonates. Furthermore, average surface CO₂ concertation available for eCO₂RR is much higher at 40 $\mu$m thick boundary layer, compared with CO₂ concertation at other boundary layer thicknesses, see Figure 8. This finding is consistent with gas flow electrolyzer using GDL, where boundary layer thickness inside the CL is only a few micrometers thick, minimizing the resistance to CO₂ mass transfer in aqueous phase [24]. Boundary layer thickness has a substantial impact on the movement of other species, such OH^-1 ions that governs pH in the vicinity of WE, and hence the selectively of eCO₂RR products [54]. Thinner boundary layer expedites the transport of OH^-1ions from the surface of CL to the bulk of an electrolyte. It is evident from Figure 8, where pH of 14.23 at boundary layer thickness of 160 $\mu$m drops to 14.08 at boundary layer thickness of 40 $\mu$m.

CL has a porous structure. Its porosity is defined as [24]:

$${\epsilon }_{CL}=1-{\displaystyle \frac{m_{loading}}{{\rho }_{CL}{L}_{CL}}}$$

(23)

$m_{loading}$ is the catalyst mass per unit area, ${\rho }_{CL}$ is the catalyst mass density, and $L_{CL}$ is the thickness of CL. Porous CL offers more active sites for eCO₂RR, than a traditional single metal solid electrode [24]. Although porous CL has more active surface area, it also minimizes catalyst loading (if the CL length is kept constant, see equation 23), which affects the TCD. Porous CL accelerates species diffusion rate; allows more CO₂ to penetrate inside CL, and the removal of product species from CL to the bulk of an electrolyte. CL has pore sizes in the range 10 – 100 nm [27]. The effect of CL porosity on TCD is plotted in Figure 9. TCD rises gradually from 50.48 mA cm^-2 to 114.32 mA cm^-2, when CL porosity is reduced from 0.8 to 0.3. The increase in TCD with reduced CL porosity is due to increase in catalyst mass loading that substantially increases active sites for CO₂ reduction. However, lower CL porosity affects CO₂ diffusion rate and increases the probability of CO₂ carbonation. Therefore, the average CO₂ concentration is lower in CL with lower porosity value, see Figure 10.

CL thickness, and CL porosity have a direct link with each other, see equation 23. Desired CL thickness with required CL porosity may be obtained by adjusting the catalyst loading. CL thickness plays a vital role in electrochemical kinetics. Although thicker CL seems to have increased catalyst loading, which could have a substantial impact on TCD; ironically, most part of the CL is rendered inoperable due to limited supply of CO₂, and fast CO₂ reduction kinetics. Moreover, thicker CL has longer tortuous diffusion path, reducing species diffusion rate inside CL. Figure 11 shows that TCD increases by reducing CL thickness. TCD with 20 $\mu$m thick CL is 71.13 mA cm^-2, and TCD with 3 $\mu$m thick CL is 93.28 mA cm^-2, approximately. A total increase of 24 mA cm^-2 is observed just by varying the CL thickness, keeping all other parameters constant. This increase in TCD is attributed to improvement in species diffusion rate. Additionally, average surface CO₂ concentration increases with decreasing CL thickness, due to shorter diffusion path inside CL, see Figure 12. Local distribution of CO₂ in computational domain (boundary layer + CL) shows that CO₂ concertation decreases gradually from the outer periphery of boundary layer to the interface of boundary layer and CL, but a sharp decline in CO₂ concertation profile is seen when CO₂ enters into the CL domain, see Figure 13. Local CO₂ concentration profile in CL shows that most CO₂ gets consumed before even reaching the middle of the CL. Because in boundary layer CO₂ undergoes a reversible chemical reaction, whereas in the CL, CO₂ is not only chemically converted but also consumed electrochemically.

Nature of an electrolyte has a significant impact on CO₂ reduction kinetics. KOH is widely recognized as one of the most effective electrolytes for CO₂ reduction, owing to its high ionic conductivity and alkaline properties, which help suppress HER [12,14,55,56,57,58]. Potassium bicarbonate (KHCO₃), a buffer solution, is known for providing near neutral pH conditions for CO₂ reduction [24,25]. To analyze the effect of nature of an electrolyte on CO₂ reduction kinetics, 1 M KHCO₃ electrolyte is simulated, and compared with 1 M KOH simulation results, see Figure 14. Two additional chemical reactions have been included in the model for 1 M KHCO₃ electrolyte, see Table 7. Kinetic rate parameters for these reactions are reported in Table 6.

There is a significant difference between the TCD obtained using KHCO₃ electrolyte and KOH electrolyte. The maximum TCD with KHCO₃ is 101.75 mA cm^-2 at an applied cathodic potential of -1.1 V vs RHE, while TCD with KOH at -1.1 V vs RHE is 116.94 mA cm^-2. This difference is due to the intrinsic nature of the electrolytic solutions. KHCO₃ is a buffer solution, therefore, to maintain the buffer it consumes more CO₂. A batch cell, already constrained by limited CO₂ supply, experiences a CO₂ shortage in the CL domain, leading to a lower TCD. A comparison of average surface CO₂ concentration for both electrolytes is shown in the Figure 15. Although the CO₂ concentration in both cases exhibits the same trend with increasing applied potential, there is an order of magnitude difference in CO₂ concentration. PCDs of CO₂ reduction products using 1 M KHCO₃ electrolyte are shown in Figure 16. It shows that electrolyte could alter the selectivity of eCO₂RR products.

The primary purpose of an electrolyte is to enhance conductivity of an aqueous phase. An electrolyte with higher molarity provides a more conductive medium by minimizing charge transfer resistance between the electrodes [45] and by reducing the negative overpotential of CO₂ reduction products [10]. Effect of KOH electrolyte strength on TCD at -1.0 V vs RHE is shown in Figure 17. TCD increases with rising electrolyte molarity but begins to plateau at significantly higher molarity levels. TCD increases as the enhanced medium conductivity facilitates the movement of charged species, thereby accelerating electrochemical kinetics. The bending of the slope at higher molarity values is attributed to the limited availability of CO₂ for eCO₂RR, see Figure 18. On one hand, KOH with higher molarity results in increased TCD; however, the higher concentration of OH⁻ ions consume more CO₂, leading to a reduction in carbon efficiency.

TCD comparison of SE, GDLE, and GCE is shown in Figure 19. There is a stark difference between TCD of SE, and the electrodes with porous structure (GDLE, GCE). TCD of SE increases with the increase in applied cathodic potential, however, the increase is not significant, in comparison with the porous electrodes. At -0.7 V vs RHE, TCD of SE is 73.94 mA cm^-2, and it rises to 77.34 mA cm^-2, at -1.1 V vs RHE. There is barely an increase of 3.39 mA cm^-2 in TCD. Furthermore, the slope of TCD curve for SE decreases with the increase in applied cathodic potential, because the availability of active surface area for CO₂ reduction becomes a limiting factor. Average surface CO₂ concentration for SE is relatively much higher compared with GDLE, and GCE, due to the limited active surface area for eCO₂RR, see Figure 20. TCD of GDLE traces the GCE curve, however, GDLE has a slightly higher TCD, as it facilitates more CO₂ to diffuse into CL. Also, average surface CO₂ concentration is higher for GDLE than GCE.

4. Conclusion

In conclusion, we have developed a numerical model to study CO₂ reduction reaction in a batch-cell with three different WEs (SE, GCE, GDLE). Numerical model is simulated using COMSOL Multiphysics software. The key conclusions derived from our numerical simulation are as follows.

• Thickness of boundary layer has an inverse relation with the TCD. Thicker boundary layer offers more resistance to CO₂ diffusion, while thinner boundary layer improves the mass transfer of the species. Magnetic stirrer speed could be manipulated to control the boundary layer thickness in a batch-cell.
• 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 current density, compared with a relatively thicker CL. Batch-cell suffers from poor CO₂ diffusion rate, and limited amount of CO₂ is available for 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 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. 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 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-perform SE in terms of current density. In particular, the porous structure of GDLE facilitates more CO₂ into the CL that increases the rate of eCO₂RR.