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

Abstract

Batch cells are pivotal in advancing the foundational research of CO₂ reduction by providing precise control over reaction conditions to study catalyst behavior and reaction mechanisms, generating insights that drive the development of scalable systems like flow reactors and ultimately supporting sustainability through the industrial adoption of carbon-neutral technologies. Therefore, a one-dimensional numerical model is developed to study electrochemical CO₂ reduction reactions in a batch cell with three different working electrode configurations: solid electrode, glassy carbon electrode, and gas-diffusion-layer electrode. The experimental results of two Cu-based catalysts are used to obtain electrochemical kinetic parameters and to validate the numerical model. The simulation results demonstrate that both gas-diffusion-layer electrodes and glassy carbon electrodes with porous catalyst layers have superior performance over solid electrodes in terms of total current density. Furthermore, we studied the impact of the key parameters of batch cells with glassy carbon electrodes, such as boundary-layer thickness, catalyst-layer thickness, catalyst-layer porosity, electrolyte nature, and the strength of an electrolyte relative to the total current density at a fixed applied cathodic potential of −1.0 V vs. RHE.

1. Introduction

The widespread use of fossil fuels in recent decades has significantly increased greenhouse gas emissions, with atmospheric CO₂ levels rising sharply from 280 ppm to 410 ppm since the industrial era [1,2,3]. The electrochemical CO₂ reduction represents a transformative approach to address pressing global sustainability challenges and aligns with global efforts to achieve carbon neutrality by closing the carbon loop and the transition to renewable energy systems [4]. When powered by renewable energy sources, such as wind or solar, these processes not only lower atmospheric CO₂ levels but also contribute to the production of sustainable fuels, and chemicals offer a pathway to recycle carbon emissions, reduce reliance on fossil-based chemicals, and promote a more sustainable society [5]. Batch cells play a critical role in the foundational research of CO₂ reduction, offering precise control over reaction conditions to study catalyst behavior, reaction kinetics, and product selectivity. Batch-cell studies are indispensable in advancing this field, enabling the discovery and optimization of catalysts, and insights gained from batch-cell studies inform the development of scalable systems, such as flow reactors [6,7,8]. Thus, CO₂ reduction in batch cells serves as a decisive step toward realizing a sustainable and circular carbon economy.

The 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. [9] by presenting a mathematical model for the estimation of local CO₂ concentrations and local pH in a boundary layer over a planar electrode at a certain current density. Hashiba et al. [10] presented a 1-D model for aqueous systems, assuming 100% CH₄ faradaic efficiency, with a 100 $\mathrm{\mu }$m thick hydrodynamic boundary layer near a working electrode (WE). This model was used to demonstrate the effect of electrolyte buffer capacities on CO₂ mass transport and pH profiles within the boundary layer. Corpus et al. [11] coupled a 1-D continuum planar electrode model with covariance matrix adaptation evolution strategies to identify the kinetic rate parameters in the Tafel equation without mass transfer limitation effects. The traditional method for obtaining kinetic rate parameters is accompanied by unquantifiable errors that lead to the overprediction of current density. Bui et al. [12] developed a 1-D transient-boundary-layer model with copper (Cu) planar electrodes to study pulsed CO₂ electrolysis. This work established that pulsing dynamically changes pH and CO₂ concentrations near solid Cu electrodes, which favors C₂₊ products. Qui et al. [13] 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₂ concentrations and their effect on product selectivity, mainly alcohols. Other than batch cells and H-cells, numerical models for CO₂ reduction in flow electrolyzers have been reported in the literature as well. Whipple et al. [14] developed a 2-D numerical model for microfluidic electrolyzers to evaluate catalysts 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 [15] that represents the cathodic section of flow electrolyzers, with silver (Ag) catalysts to demonstrate the superiority of GDEs over planar electrodes; and a zero-gap MEA model [16], with Ag catalysts to study gas-fed (full-MEA) and liquid-fed (exchange-MEA) configurations. This model was extended to include Cu kinetics [17]. Kas et al. [18] developed a 2-D numerical model for GDEs to emphasize the concentration gradients in the flow channels. Ehlinger et al. [19] employed a 1-D planar electrode model to determine the kinetics of hydrogen formations from bicarbonates. The kinetic parameters obtained from the planar electrode model were then implemented in a 1-D MEA model to investigate flooding and to conduct a sensitivity analysis of MEA designs and operating parameters. In addition, in numerous studies, numerical modeling is used to support experimental findings such as carbonation, water management, local pH and species concentration, and electrolyzer optimization [20,21,22,23,24,25,26,27,28,29].

The 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 the literature on the numerical modeling of batch cells with WEs other than planar electrodes, whereas most researched electrocatalysts are available in nano-particle forms [4,30,31] that require some support in order to be tested in an electrolyzer. Gas-diffusion-layer electrodes (GDLEs) and glassy carbon electrodes (GCEs) are employed to support electrocatalysts that are available in nano-particle form [32]. The aforementioned references in the literature are about planar electrodes (solid electrodes) in batch cells, H-cells, and flow-type H-cells; and GDEs in H-cells with flow configurations, in microfluidic cells, in flow electrolyzers with catholytes, 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 cells with three different Wes, including solid electrodes (SEs), GCEs, and GDLEs [33,34]. The numerical model is validated with the experimental results of CuAg_0.5Ce_0.2 catalysts developed indigenously in our lab in a batch cell with GCEs. The numerical model is further tested with the published experimental results of Ebaid et al. [35] for CO₂ reduction over Cu foil in a gas-tight electrolyzer with anion exchange membranes separating cathodes and anodes. The experimental setup details of the former catalyst including catalyst synthesis, batch-cell operation, and result analysis are reported in the Supplementary Materials. The numerical model is then used to investigate total current density (TCD) behavior with respect to variations in batch-cell parameters, including boundary-layer thickness, catalyst-layer (CL) porosity, CL thickness, the nature of an electrolyte, and the electrolyte’s strength. This study provides information for experimentalists regarding the selection of WEs and a critical analysis of the key parameters associated with batch cells that could be adjusted for the desired performance.

Although the numerical model demonstrates complex physical and chemical processes inside the boundary layer and the porous CL (GCE and GDLE) and it shows good agreement with experimental findings, some assumptions have been adopted in the model that should be replaced in future works with relevant physics to capture even more realistic results. In batch-cell experiments, a 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 the bulk of an electrolyte and boundary layer. The boundary-layer thickness is a strong function of stirring speeds. Higher stirring speeds result in a thinner boundary layer, and the vice versa case is possible, but in the model, the boundary-layer thickness is adopted from the literature [6,12,15,18,36], and the impact of various boundary-layer thicknesses has been discussed in detail. Cu-based electrocatalysts are known for producing a range of gaseous and liquid CO₂ reduction products. During batch-cell operations, there is a high probability that some CO₂ reduction products interact with each other and influence pH, CO₂ solubility, and species transport. To avoid 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 the applied potential and the heat from 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 cells is the WE of GCE, where CO₂ is reduced. GCE has a solid flat surface coated with CL. The boundary layer is formed on the surface of the catalyst that governs the mass transfer of species between the bulk of an electrolyte and WE, while the bulk of an electrolyte achieves uniformity due to constant stirring. Therefore, the computational domain for GCE has two important regions: the porous CL and the boundary layer. A schematic of the one-dimensional (1-D) computational domain is shown in Figure 1 with all the necessary boundary conditions.

A comprehensive numerical model is developed to deconvolute physio-electrochemical phenomena 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 were modeled: dissolved CO₂, hydroxyl ions (OH⁻¹), bicarbonate ions (HCO₃⁻¹), carbonate ions (CO₃²⁻), protons (H⁺), and potassium ions (K⁺). The 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

Table 1. Electrochemical CO₂ reduction reactions for CuAg_0.5Ce_0.2 catalysts.

Species	Reactions
H2:	$2{\mathrm{H}}_2{\mathrm{O}}_{\left(\mathrm{l}\right)}+2{\mathrm{e}}^{-}\to {{\mathrm{H}}_2}_{\left(\mathrm{g}\right)}+2\mathrm{O}{\mathrm{H}}^{-}$
CO:	$C{O_2}_{\left(l\right)}+H_2{O}_{\left(l\right)}+2e^{-}\to C{O}_{\left(g\right)}+2O{H}^{-}$
CH4:	$C{O_2}_{\left(l\right)}+6H_2{O}_{\left(l\right)}+8e^{-}\to C{H_4}_{\left(g\right)}+8O{H}^{-}$
C2H4:	$2C{O_2}_{\left(l\right)}+8H_2{O}_{\left(l\right)}+12e^{-}\to C_2{H_4}_{\left(g\right)}+12O{H}^{-}$
EtOH:	$2C{O_2}_{\left(l\right)}+9H_2{O}_{\left(l\right)}+12e^{-}\to C_2{H}_{5}O{H}_{\left(l\right)}+12O{H}^{-}$
PrOH:	$3C{O_2}_{\left(l\right)}+13H_{2}O+18e^{-}\to C_3{H}_{7}O{H}_{\left(l\right)}+18O{H}^{-}$

[30]. It is assumed that gaseous products (H₂, CO, CH₄, and 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 a ubiquitous species in the computational domain; therefore, it is not included in modeling. Electrolyte is a 1 M KOH solution. Besides the electrochemical conversion of CO₂ into useful products, a considerable amount of CO₂ reacts reversibly with the hydroxyl ions (OH⁻¹) to produce bicarbonate (HCO₃⁻¹) and carbonate (CO₃²⁻) ions. These reactions govern the local pH in the computational domain. Due to the high pH (~1 M KOH pH = 14) of an electrolyte, only the reactions in

Table 2. Chemical reaction of CO₂ with electrolytes in alkaline media.

R1:	$C{O_2}_{\left(l\right)}+O{H}^{-}\rightleftharpoons HC{O}_3^{-}$	K1 (1/M)
R2:	$HC{O}_3^{-}+O{H}^{-}\rightleftharpoons C{O}_3^{2-}+H_{2}O$	K2 (1/M)
Rw:	$H_{2}O\rightleftharpoons H^{+}+O{H}^{-}$	Kw (M2)

have been considered in the model [37]. The precipitation of bicarbonate/carbonate salts is ignored in the numerical model. Although salt build-up has a substantial deteriorating effect on the CL [38,39], it also increases the required cell potential if the experiment is carried out for an extended period [40], but in the case of batch-cell modeling, this phenomenon can be ignored as a fresh electrolyte solution is used for each batch-cell experiment [38].

Species balance is used to model the 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 species j in the computational domain, while source terms on the right side represent the consumption and generation of species due to chemical and electrochemical reactions, respectively. The Nernst–Planck equation is invoked for the molar flux of species in the computational domain [34].

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

(2)

The 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 the Bruggeman correlation [41]. 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 the electrolyte potential (${\varphi }_l$) to close the degree of freedom; therefore, an electroneutrality equation is included in the model.

$$\sum _{j=1-6}{z}_j{c}_j=0$$

(5)

The Nernst-Planck equation assumes dilute-solution theory [15,34]. The modeled domain is not dilute under specified conditions. The concentrated solution theory requires additional diffusion coefficients that 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 species balance is defined as follows:

$$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 the molar weight of species. $a_s$ is the effective surface area available for CO₂ reduction in porous CL domains, $n_k$ represents the electron transferred in each reaction $k$, $i_k$ is the partial current density (PCD), and F is Faraday’s constant. The electrode’s kinetics are modeled using the concentration-dependent Tafel equation [42], which describes the electrochemical conversion of dissolved CO₂. A detailed derivation of the Tafel equation is provided in the Supplementary Materials. The final form of the Tafel kinetic equation used in the model is as follows:

$$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 the Tafel kinetic equation is adopted from Bui et al. [12]. The Tafel equation has several unknown parameters that are either the direct inputs to a simulation or obtained through data fitting to experimental findings. ${\varphi }_s$ is the applied cell potential, while ${\varphi }_l$ is evaluated by solving the Nernst–Planck equation and electroneutrality equation (see Equations (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 the experimental data fitting of CO₂ reduction at various applied potentials [16] or through rigorous microkinetic modeling [36]. According to density functional theory calculations, ${\gamma }_{C{O}_2}$ is the number of elementary steps in CO₂ reduction toward certain product species, which is the most simplified way to obtain this parameter. ${\gamma }_{pH}$ is the pH rate ordering parameter. It is estimated by fitting the experimental data of the 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 reference [12], and their values are listed in Table 4. $c_{C{O}_2}$ and $pH$ are local CO₂ concentrations and local pH, respectively. Obtaining these values experimentally is quite challenging; hence, simulations are conducted, similarly to Gupta et al. [9], 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 the computational domain is modeled using the following equations.

$$\nabla ·i=0$$

(10)

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

(11)

$$S=-a_s\sum _{k=1-6}{i}_k$$

(12)

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

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

(13)

${\sigma }_m$ represents the conductivity of solid catalysts and liquid electrolytes. In porous CL, effective conductivity is used according to the Bruggeman correlation.

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

(14)

The 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 follows [43,44]:

$$\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 the case of an electrolyte with a high concentration of ions, the solubility of CO₂ is reduced. Reduced CO₂ solubility in an electrolytic solution is evaluated using Sechenov’s constants [45].

$${\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 constant values are provided in

Table 3. List of Sechenov’s constants [45].

Constants	Values
$h_{G,0}$	−0.0172
$h_T$	−0.000338
$h_K$	0.0922
$h_{OH}$	0.0839
$h_{HC{O}_3}$	0.0967
$h_{C{O}_3}$	0.1423

.

$$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 the 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 the intricate network of pores, the diffusion of CO₂ inside CL is reduced. The nature of eCO₂RR products depends on the catalyst. The thickness/length of the boundary layer is about 100 $\mathrm{\mu }\mathrm{m}$ [12]. However, this value changes with the stirring speed (in RPM). The potential (${\varphi }_s=E_{WE}$) is applied at x = 0, whereas it is assumed that the reference electrode (RE) is located at the interface of the boundary layer and bulk electrolyte (x = L_CD), where the electrolyte potential is set to zero (${\varphi }_l=0$). Moreover, steady-state and isothermal conditions are assumed.

2.1. Kinetic Parameters from Experimental Data

In this section, the Tafel equation parameters for each CO₂ reduction product are obtained from experimental data. The experimental data are provided in Table S1 (refer to Supplementary Materials). The 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}$, and ${\alpha }_{c,k}$. The parameter values for ${\gamma }_{C{O}_2,k}$ and ${\gamma }_{pH,k,SHE}$ are provided in

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

Product Species	${\mathit{\gamma }}_{\mathit{C}{\mathit{O}}_{\mathbf{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

. Parameters $c_{C{O}_2}$ and $pH$ are listed in

Table 5. Pre-simulation estimates for local CO₂ concentrations and pH values using the experimental data of CuAg_0.5Ce_0.2 catalysts.

TCD	CO2	pH
(mA cm−2)	(mM)
−47	15.747	11.072
−64	12.150	11.128
−84	8.0113	11.219
−105	6.1246	11.28
−131	4.2034	11.368

. $i_{0,k}$, and ${\alpha }_{c,k}$ values are reported in

Table 6. Batch-cell model parameters.

Parameter	Value	Unit	Reference
Design Parameters
LCL	15	$\mathrm{\mu }\mathrm{m}$	Experimental input
LBL	100	$\mathrm{\mu }\mathrm{m}$	Fitted
LGDL	325	$\mathrm{\mu }\mathrm{m}$	Experimental input
T	298	K	This study
EWE	−0.7–−1.1	V vs. RHE
Electrochemical Kinetic Rate Parameters for CuAg0.5Ce0.2 Catalysts
$i_{0,h2}$	97.3	mA cm−2	This Study
$i_{0,CO}$	8.89 × 108	mA cm−2
$i_{0,C{H}_4}$	9.50 × 105	mA cm−2
$i_{0,C_2{H}_4}$	2.8	mA cm−2
$i_{0,EtOH}$	4.10 × 10−2	mA cm−2
$i_{0,PrOH}$	9.50 × 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,38]
$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 × 10−3	m3 mol−1s−1	[18]
$k_{2,f}$	1.0 × 10−8	m3 mol−1s−1	[18]
$k_{3,f}$	0.04	s−1	[15]
$k_{4,f}$	56.281	s−1	[15]
$k_{w,f}$	0.001	mol L−1s−1	[17]
$k_{1,b}$	1.34 × 10−4	s−1	[18]
$k_{2,b}$	2.15 × 10−4	s−1	[18]
$k_{3,b}$	9.37 × 104	L mol−1s−1	[15]
$k_{4,b}$	1.23 × 1012	L mol−1s−1	[15]
$k_{w,b}$	1.0 × 1011	L mol−1s−1	[17]
Porous Media Properties
${\epsilon }_{GDL}$	0.72		This Study
${\epsilon }_{CL}$	0.3		This Study
${\sigma }_{GDL}$	220	S m−1	[15]
${\sigma }_{CL}$	100	S m−1	[15]
$a_s$	1.0 × 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	[46]
$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	[46]
$D_{HC{O}_3^{-}}$	$7.016\times {10}^{-9}{\left({\displaystyle \frac{T}{204.0282}}-1\right)}^{2.3942}$		m2 s−1	[47]
$D_{C{O}_3^{=}}$	$5.447\times {10}^{-9}{\left({\displaystyle \frac{T}{210.2646}}-1\right)}^{2.1929}$		m2 s−1	[47]
$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	[48]
$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	[48]

.

The local CO₂ concentration and pH values are estimated by simulating experimental data. The experimental TCD is set as a boundary condition. The mass flux boundary condition, estimated from experimental PCD, is used for dissolved CO₂ and hydroxyl ion (OH⁻¹) species balance. In this way, the local average CO₂ concentration and local average pH values are estimated in the computational domain. These values are provided in Table 5.

The remaining two parameters, $i_{0,k}$ and ${\alpha }_{c,k}$, are obtained through comparing with experimental data. The 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)

The Tafel equation is simplified as follows:

$$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 Equation (22) are denoted as ${\alpha }_{c,k}$ and $i_{0,k}$, respectively. The normalized current density plots as a function of the overpotential for each experimentally recorded CO₂ reduction product are shown in Figure S3.

2.2. Numerical Model for SE and GDLE

In this section, the numerical models for the other two WEs are presented: SE and GDLE. In practice, SE is made up of a single metal, such as Cu, but in this study, it is assumed that SE is made up of the catalyst CuAg_0.5Ce_0.2. GDLE is a GDL that contains both macropores (carbon fibers) and micropores (carbon powder). A layer of catalysts is coated on top of the micropore 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 solely comprises the boundary layer. The computational domain contains six species: dissolved CO₂, OH⁻, HCO₃⁻, CO₃⁼, H⁺, and K⁺. Species flux is solved using Equation (1). The concentration of these species is set to the bulk electrolyte values at (x = L_CD). CO₂ diffuses from the interface of the boundary layer and bulk electrolyte (x = L_CD) to reach the surface of the electrode. Electrochemical reactions happen on the surface of the electrode (x = 0). Therefore, the 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). The potential (${\varphi }_s$) is applied at the surface of the electrode, while the electrolyte potential (${\varphi }_l$) is set to 0 at x = L_CD, assuming that RE is placed at this location. In the case of SE, the length of the computational domain is similar to the thickness of the boundary layer, and the length of the computational domain is 100 $\mathrm{\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). The GDLE computational domain has 4 regions: 2 boundary layers, a porous GDL, and a porous CL. The computational domain contains dissolved CO₂, OH⁻, HCO₃⁻, CO₃⁼, H⁺, and K⁺, and species transport is modeled using Equation (1). The 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. The GDL serves two important functions: It supports the CL and acts as a medium for electronic conductors. The potential (${\varphi }_s$) is applied at the interface of the 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

The numerical model developed in this study is simulated using the COMSOL Multiphysics v6.1 environment. COMSOL Multiphysics uses finite element techniques to discretize the differential equations. A MUMPS direct general solver is used to solve a set of algebraic equations (Ax = B). The total number of nodes in the computational domain is 200; the number of nodes in the CL domain and the boundary-layer domain is set to 150, and 50, respectively. The information required to simulate the model is provided in Table 6.

A grid-independent study is performed to check the model’s consistency and robustness. The CL domain is more sensitive to the species concentration gradient than the boundary-layer domain, as eCO₂RR happens inside the pores of the CL; therefore, a smaller element size (0.30 $\mathrm{\mu }$m) is used in the CL domain, comparing the boundary-layer domain with the relatively bigger element size (0.67 $\mathrm{\mu }$m). A 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 the number of nodes (decreasing element size) in the computational domain. The grid-independent solution is shown in Figure 4.

2.4. Numerical Model Validation

The numerical model is validated by comparing the simulated TCD result with experimental data (see Figure 5). The numerical model shows good agreement with the experiment. However, there is a slight difference between the simulated TCD slope and the experimental TCD slope. The numerical model predicts TCD with a positive slope; however, the slope of the experimental data shows a decreasing trend with an increase in applied cathodic potential. In experiments, electrolytes saturated with CO₂ are used, and the batch cell is sealed for the duration of the experiment. The 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 to 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 experiments where mass transfer resistance increases with decreasing CO₂ concentrations in the bulk electrolyte. Therefore, a decreasing slope trend for the experimental current density is due to the scarce supply of CO₂, which becomes significant at higher applied potentials. Furthermore, the use of raw techniques for kinetic parameter estimation and the incapability of the 1-D model to fully describe the underlying physics could also affect the numerical simulation’s results. Despite the simplification of the numerical model, it is able to predict TCD with quite good accuracy.

A comparison of the simulated PCD of eCO₂RR products with experimental data is shown in Figure S4. Except for ethylene, all individual current densities show a good match with the experiment. The experimental data for ethylene show irregular behavior, especially the data point at −0.7 V vs. RHE and −1.0 V vs. RHE, which significantly depart from the normal trend.

The numerical model developed in this study is also validated with the experimental data of Ebaid et al. [35]. A Cu foil and graphite rod were used as a WE and counter electrode (CE), respectively. The electrolyte was 0.1 M CsHCO₃ saturated with CO₂. The validation results for TCD and PCD are shown in Figure 6 and Figure S5, respectively.

3. Results and Discussion

In this section, a numerical model is used to study the effect of the critical parameters of batch cells with GCEs using the kinetic parameters of CuAg_0.5Ce_0.2 catalysts. The thickness of the boundary layer is a crucial factor that governs the diffusion rate of species. Thicker boundary layers offer significant resistance to CO₂ diffusion. They do not only increase the diffusion time of CO₂ to reach the catalyst’s surface. The probability of CO₂ consumption via side reactions also increases, even further reducing the 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 thicknesses on TCD at −1.0 V vs. RHE is shown in Figure 7. The TCD is obtained for four different boundary-layer thicknesses of equal intervals: 40 $\mathrm{\mu }$m, 80 $\mathrm{\mu }$m, 120 $\mathrm{\mu }$m, and 160 $\mathrm{\mu }$m. Figure 7 shows an inverse relationship between the TCD and boundary-layer thickness, and the trend is not linear. The substantial improvement in the TCD at a boundary-layer thickness of 40 $\mathrm{\mu }$m, compared with 160 $\mathrm{\mu }$m, is due to the lower resistance to CO₂ mass transfers and the decline in the chemical conversion of CO₂ into bi-carbonates and carbonates. Furthermore, the average surface CO₂ concentration available for eCO₂RR is much higher at the 40 $\mathrm{\mu }$m thick boundary layer compared with the CO₂ concentration at other boundary-layer thicknesses (see Figure 8). This finding is consistent with gas flow electrolyzers that use the GDL, where the boundary-layer thickness inside the CL is only a few micrometers thick, minimizing the resistance to CO₂ mass transfers in the aqueous phase [15]. The boundary-layer thickness has a substantial impact on the movement of other species, such as the OH⁻¹ ions that govern the pH in the vicinity of WE, hence the selectivity of eCO₂RR products [49]. A thinner boundary layer expedites the transport of OH⁻¹ions from the surface of the CL to the bulk of an electrolyte. This is evident in Figure 8, where the pH of 14.23 at a boundary-layer thickness of 160 $\mathrm{\mu }$m drops to 14.08 at a boundary-layer thickness of 40 $\mathrm{\mu }$m.

The CL has a porous structure. Its porosity is defined as follows [15]:

$${\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 the CL. The porous CL offers more active sites for eCO₂RR than a traditional single-metal solid electrode [15]. Although the porous CL has more active surface areas, it also minimizes catalyst loading (if the CL’s length is kept constant, see Equation (23)), which affects the TCD. The porous CL accelerates species diffusion rates, allows more CO₂ to penetrate inside the CL, and removes the product species from the CL to the bulk of an electrolyte. The CL has pore sizes within the range of 10–100 nm [18]. The effect of CL porosity on TCD is plotted in Figure 9. TCD rises gradually from 50.48 mA cm⁻² to 114.32 mA cm⁻², when the CL’s porosity is reduced from 0.8 to 0.3. The increase in the TCD with reduced CL porosity is due to an increase in catalyst mass loading that substantially increases the active sites for CO₂ reduction. However, a lower CL porosity affects the CO₂ diffusion rate and increases the probability of CO₂ carbonation. Therefore, the average CO₂ concentration is lower in CL, with lower porosity values (see Figure 10).

CL thickness and CL porosity have a direct link with each other (see Equation (23)). The desired CL thickness with respect to the required CL porosity may be obtained by adjusting the catalyst’s loading. The CL’s thickness plays a vital role in electrochemical kinetics. However, a thicker CL seems to have increased catalyst loading, which could have a substantial impact on TCD; ironically, most parts of the CL are rendered inoperable due to the limited supply of CO₂ and fast CO₂ reduction kinetics. Moreover, a thicker CL has longer tortuous diffusion paths, reducing the species diffusion rate inside the CL. Figure 11 shows that the TCD increases by reducing CL’s thickness. The TCD with a 20 $\mathrm{\mu }$m thick CL is 71.13 mA cm⁻², and the TCD with a 3 $\mathrm{\mu }$m thick CL is 93.28 mA cm⁻² approximately. A total increase of 24 mA cm⁻² is observed just by varying the CL’s thickness, keeping all other parameters constant. This increase in the TCD is attributed to improvements in the species diffusion rate. Additionally, the average surface CO₂ concentration increases with decreasing CL thicknesses due to the shorter diffusion path inside the CL (see Figure 12). The local distribution of CO₂ in the computational domain (boundary layer + CL) shows that CO₂ concentrations decrease gradually from the outer periphery of the boundary layer relative to the interface of the boundary layer and CL, but a sharp decline in the CO₂ concentration profile is seen when CO₂ enters into the CL domain (see Figure 13). Local CO₂ concentration profiles in the CL show that most CO₂ become consumed before even reaching the middle of the CL. This is because CO₂ undergoes a reversible chemical reaction in the boundary layer, whereas in the CL, CO₂ is not only chemically converted but also consumed electrochemically.

The 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 due to its high ionic conductivity and alkaline properties, which help suppress HER [22,38,50,51,52,53]. Potassium bicarbonate (KHCO₃), a buffer solution, is known for providing near-neutral pH conditions for CO₂ reduction [24,25]. To analyze the effect of the 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 the 1 M KHCO₃ electrolyte (see

Table 7. Additional reversible chemical reactions.

R3:	$C{O_2}_{\left(l\right)}+H_{2}O\rightleftharpoons H^{+}+HC{O}_3^{-}$	K3 (M)
R4:	$HC{O}_3^{-}\rightleftharpoons H^{+}+C{O}_3^{2-}$	K4 (M)

). The kinetic rate parameters for these reactions are reported in Table 6.

There is a significant difference between the TCD obtained using KHCO₃ electrolytes and KOH electrolytes. The maximum TCD with KHCO₃ is 101.75 mA cm⁻² 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⁻². 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 the limited CO₂ supply, experiences a CO₂ shortage in the CL domain, leading to a lower TCD. A comparison of average surface CO₂ concentrations for both electrolytes is shown in 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₂ concentrations. The PCDs of CO₂ reduction products using 1 M KHCO₃ electrolyte are shown in Figure 16. They show that electrolytes could alter the selectivity of eCO₂RR products.

The primary purpose of an electrolyte is to enhance the conductivity of an aqueous phase. An electrolyte with higher molarity provides a more conductive medium by minimizing the charge transfer resistance between the electrodes [37] and reducing the negative overpotential of CO₂ reduction products [40]. The effect of the KOH electrolyte’s strength on the TCD at −1.0 V vs. RHE is shown in Figure 17. The 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 the one hand, the KOH with higher molarity results in an increased TCD; however, the higher concentration of OH⁻ ions consumes more CO₂, leading to a reduction in carbon efficiency.

A TCD comparison of SE, GDLE, and GCE is shown in Figure 19. There is a stark difference between the TCD of SE and the electrodes with porous structures (GDLE and GCE). The TCD of SE increases with an increase in the applied cathodic potential; however, the increase is not significant in comparison with porous electrodes. At −0.7 V vs. RHE, the TCD of SE is 73.94 mA cm⁻², and it rises to 77.34 mA cm⁻² at −1.1 V vs. RHE. There is barely an increase of 3.39 mA cm⁻² in TCD. Furthermore, the slope of the TCD curve for SE decreases with an increase in the applied cathodic potential because the availability of the active surface area for CO₂ reduction becomes a limiting factor. The 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). The TCD of GDLE traces the GCE curve; however, the GDLE has a slightly higher TCD as it facilitates more CO₂ to diffuse into the CL. Also, the average surface CO₂ concentration is higher for the GDLE than the GCE.

4. Conclusions

In conclusion, we developed a numerical model to study CO₂ reduction reactions in a batch cell with three different WEs (SE, GCE, and GDLE). A numerical model was simulated using COMSOL Multiphysics v6.1 software. The key conclusions derived from our numerical simulation are as follows:

• The thickness of the boundary layer has an inverse relation with the TCD. A thicker boundary layer offers more resistance to CO₂ diffusion, while a thinner boundary layer improves the mass transfer of the species. The magnetic stirrer’s speed could be manipulated to control the boundary layer’s thickness in a batch cell.
• CL porosity also demonstrates an inverse relation with the TCD; the higher the CL’s porosity, the lower the TCD, 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 a good TCD.
• The CL’s 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. The batch cell suffers from poor CO₂ diffusion rates, and a 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 the catalyst, improves the 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.
• The electrolyte’s nature impacts the CO₂ reduction reaction. Moreover, 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.
• The KOH electrolyte with higher molarity increases the TCD by minimizing the charge transfer resistance between electrodes. Nevertheless, the 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 electrolytes for conductivity improvement.
• Among the three electrodes, the GDLE and GCE outperform the SE in terms of current density. In particular, the porous structure of the GDLE facilitates more CO₂ into the CL, which increases the rate of eCO₂RR.