Sensitivity Analysis of One-Dimensional Multiphysics Simulation of CO2 Electrolysis Cell

Electrochemical (EC) carbon dioxide (CO2) reduction, where CO2 is converted to value-added products such as fuel precursors, plays a key role in helping the world’s energy system reach net-zero carbon emissions. Simulations of EC cells provide valuable insight into their operation since detailed experimental results on short length and time scales are difficult to obtain. In this work, we construct a 1D simulation of a membrane-electrode-assembly EC cell for CO2 reduction, using a porous silver gas diffusion cathode. We run the simulation under different electrolyte conditions, showing how the cell performance is affected. We then perform a sensitivity analysis of all input parameters to the simulation, which has not been presented before in the literature. We show that the CO partial current density (iCO) is significantly affected by each input parameter of the simulation. iCO is most sensitive to EC kinetic parameters (i0/α) of all EC reactions, with a 1% change in α resulting in up to 6% change in iCO. Since there is uncertainty associated with the value of each input parameter, this indicates that infidelity between experiment and simulation is likely, and thus, caution should be practiced when comparing experimental results to simulation results. Further, we show that the large range of conditions simulated in literature helps to explain the large variance in reported values of i0 and α. The results of this paper demonstrate the potential of sensitivity analysis methods to quickly optimize aspects of cell performance (CO2 utilization, Faradaic efficiency, etc.).


■ INTRODUCTION
Context.The field of electrochemical (EC) CO 2 reduction is concerned with converting CO 2 to value-added products by means of EC reactions.The field is undergoing a surge in popularity as increasing focus is shifting to net-zero-emission energy systems for the world in light of the climate crisis.EC CO 2 reduction creates a source of renewable hydrocarbons, eliminating the need for fossil resources, if the value-added products are fuels or fuel precursors.It incentivizes carbon capture by turning captured CO 2 to value-added products, and it is also a way of storing excess renewable energy in the form of hydrocarbons for later use. 1 The operation of a CO 2 electrolysis cell is complex, as gas flow, liquid flow, EC reactions, solution-phase chemical reactions, phase change of species between gas and liquid phases, electron transport, and adsorption of species onto the catalyst surface all occur simultaneously and influence each other.For this reason, it is difficult to operate these cells consistently and reliably at high current densities with high efficiencies. 2,3To aid with this problem, simulations of CO 2 electrolyzers have been increasing in number, as can be seen in the review of Bui et al. 4 Such simulations, or "digital twins" to experiment, are useful as they provide insight at length and time scales which are impossible or difficult to attain by experiment.The purpose of such simulations is to characterize the output of the cell, typically the current−voltage behavior, and the quantities of chemical species in the cell as a function of time and operating conditions.
−13 MEA cells with a GDE cathode have been shown to produce higher current densities primarily due to higher rates of mass transport to the cathode surface caused by a shorter diffusion path for CO 2 to the electrode surface. 11iterature Review.Due to constraints of time and computational power, multiphysics simulations (which are mostly in 1D) simulate some, but not all, of the physics occurring inside a real electrolysis cell.Simulations generally consider (i) transport of liquid-phase species, (ii) transport of gas-phase species, (iii) transport of electrons in solid phase,  There is a large range in the reported i 0 /α values, corresponding to the large range in simulated conditions.

The Journal of Physical Chemistry C
(iv) EC reaction kinetics at the electrode surface, and (v) solution-phase reaction kinetics in the electrolyte. 4The equations used are implemented often in the same simulation platform, COMSOL Multiphysics, although other platforms are available, e.g., ANSYS-Fluent.From these numerical frameworks, cell behavior under certain conditions and using different setups have been assessed and hypothesized.For instance, Singh et al. 6 investigate different sources of voltage loss in a photovoltaic (PV) cell and propose optimal device configurations to minimize them�suggesting the use of an AEM to minimize overall potential losses, as this allows transport of anions which are the primary charge carriers in this system.This finding is backed up by Weng et al. 11 who simulate three different cases: a full MEA (gas feed at anode and cathode) and an exchange MEA (liquid feed at anode, gas feed at cathode) using either KOH or KHCO 3 as the electrolyte.Wang et al. 14 also simulate this setup and show that the CO 2 conversion efficiency is limited to ∼50%, due to crossover of carbonate anions across the membrane from cathode to anode and subsequent conversion to CO 2 , which leaves the cell unreacted.CO 2 reduction simulations are growing in sophistication, following the trend of water reduction simulations. 46][7][8]10,11,13 Some 2D simulations of half-cells have been developed, 17−23 which give deeper insight into cell behavior as they allow us to examine perpendicular flow, as well as distribution of current and species in two dimensions. Many studis do not consider a full set of solution-phase reactions, as this presents the challenge of knowing the reaction kinetic constants as well as making model convergence more difficult.4 Simplifying assumptions have been made, for instance, Moore et al. 18 impose a concentrated KOH solution with pH above 12 and thus neglect the carbonate/bicarbonate buffer reactions (reactions R4, R5 and R7 in this work) and consider consumption of CO 2 by OH − (reaction R6) to be irreversible.The inclusion of a full solution-phase reaction system is advised, as carbonate/bicarbonate anions are identified by Rabinowitz and Kanan, 24 as well as Wang et al., 14 as limiting  2, and they are solved using COMSOL Multiphysics 6.0. Equations are described in Newman. 33 The Journal of Physical Chemistry C cell performance by reducing the carbon efficiency as will be discussed in the next section.
One key challenge in building a multiphysics simulation of a CO 2 electrolysis cell is sourcing the input parameters.Many input parameters needed for multiphysics simulations are difficult and time-consuming to measure (e.g., active surface areas of porous electrodes).For this reason, those building simulations often source parameters from work done by other groups.However, the conditions under which these parameters were originally determined do not always exactly match those of the system being simulated.This will lead to "infidelity" in the simulation�the output of the simulation will be different to that of the experiment to which it is being compared.Further, many of the input parameters commonly used may be inaccurate.For example, those who do include solution-phase reactions generally use the values presented by Bui et al., 4 or Schulz et al. 25 for their rate constants and do not account for the fact that rate constants change with conditions, e.g., under electric fields.It is common practice to round these parameters to the nearest order of magnitude from their reference source.The problem of sourcing parameters becomes more difficult since many parameters are highly variable.One example is the EC kinetic parameters, i 0 and α, used in the Butler−Volmer equations, which are used to calculate current density of an EC reaction.For every reaction, there is a large variance in the reported values for i 0 and α.Table 1 shows the reported i 0 and α values for reaction R1, CO 2 reduction to CO.These parameters are usually extracted by fitting the Butler−Volmer equations to experimental data.This is an overly simplistic method however as it does not consider effects such as mass transport, which have affected the experimental data, as pointed out by Corpus et al. 26 They suggest fitting the output of a multiphysics simulation to experimental data, thus accounting for mass transport effects and effects of solutionphase reactions.
Published multiphysics simulations are summarized in Table 1.Approximately one-third are not validated experimentally, 6,10,11,17,18,21,27 as it is difficult to obtain experimental data aside from polarization curves or Faradaic efficiencies.One-third of simulations validate their simulation using data from other published experimental setups. 8,9,12,13,15,28,29These setups, however, may not be representative of the simulated system.−32 Table 1 shows the range in reported EC kinetic parameters i 0,1 and α c,1 .Note that each case simulates a different cell setup under different operating conditions, thus corresponding to different input parameters.The choice of cell configuration, electrolyte, membrane, full or half-cell, and reactions considered are all different and thus help to explain the large discrepancy in the reported values of i 0,1 and α c,1 .Each different condition results in a different set of input parameters used by the simulation, which affects the output of the simulation and, thus, the fitting process.
In the Supporting Information, we include a further discussion of each input parameter, where they are usually sourced, and how they may lead to infidelity in the output of a simulation and thus to extraction of inaccurate kinetic parameters.
Objectives of This Work.In this work, we aim to investigate the effects of each input parameter on the current density of the CO 2 reduction reaction, i CO .To this end, we construct a multiphysics simulation of a full EC cell for CO 2 reduction, considering chemical reactions and mass transport.We consider EC reactions, reversible solution-phase reactions, and solvation of CO 2(g) to CO 2(aq) .Mass transport considered is of chemical species in the gas and solution phases and of electrons in the solid phase.We then perform sensitivity analysis of the simulation input parameters to assess how they affect the output of the simulation.This sensitivity analysis helps to explain the large variation in reported values of i 0 and α.It also suggests where more caution should be observed in sourcing input parameters to simulations of this type to minimize infidelity between experiment and simulation.Sensitivity analysis also informs future experiments by telling us what we can change in the cell to improve performance most easily.
■ METHODOLOGY Overview of Simulation.The one-dimensional multiphysics simulation presented in this work is composed of three domains as shown in Figure 1: cathode, anion-exchange membrane, and anode.The simulation accounts for EC (surface) reactions, solution-phase reactions, solvation of CO 2(g) to CO 2(aq) , transport of gaseous species (CO 2(g) , H 2(g) , CO (g) , N 2(g) ) and solution-phase species (CO 2(aq) , H 2 O (aq) , H (aq) + , OH (aq) − , HCO 3(aq) − , CO 3(aq) 2− , and K (aq) + ), and transport of electrons in the solid phase.Both electrodes are described as porous electrodes, meaning reactions and transport occur simultaneously in these domains.In the cathode domain, transport of gas-and solution-phase species occurs alongside reaction kinetics, as the cathode is a GDE.In the anode, transport of solution-phase species occurs alongside reaction kinetics.The equations in Figure 1 applied to each domain, governing transport and kinetics, are outlined below and are also outlined in the review of Bui et al. 4 The equations are populated with parameters in Table 2 and solved using COMSOL Multiphysics v6.0.
Electrochemical, Solution-Phase, and Solvation Reactions Simulated.The EC reactions occurring at the cathode are the formation of CO (g) and the unwanted formation of H 2(g) (R1) At the anode, the formation of O 2d (g) occurs due to water oxidation, which occurs via two pathways (R3a) These EC reactions are commonly used to describe CO (g) , H 2(g) , and O 2(g) formation on Ag and IrO 2 catalysts, 4 and are commonly referred to as CO 2 R, HER, and OER reactions, respectively (i.e., CO 2 reduction reaction, hydrogen evolution reaction, oxygen evolution reactions).As is the standard in multiphysics simulations, these reactions are approximated as a single step.This approximation is a source of infidelity between experiment and simulation 8 (see discussion in Results).Note that this study excludes the acidic pathway of the HER reaction [2H + + 2e − → H 2 ] due to the alkaline conditions simulated and the resulting low H + concentration.The solution-phase reactions simulated are the bicarbonate/carbonate system The Journal of Physical Chemistry C reactions (reactions R4−R7) and the dissociation of H 2 O (aq) (reaction R8), listed below: The solvation of CO 2(g) to CO 2(aq) within the cathode is also simulated: CO 2 is the only species existing in both gas and solution phases in the simulation�gas-phase water and effects, e.g., those due to humidity, are excluded since the temperature of the system is 298 K.The solution-phase reactions outlined above are those which are commonly included in multiphysics simulations of CO 2 reduction cells and whose apparent rate constants have been studied, 25 although we note that there may be reactions occurring and influencing cell performance, which are not listed here.In the Supporting Information, we include a list of such reactions.
It is suggested that the CO 2 utilization efficiency of a zerogap MEA cell with anion-exchange membranes is limited to ∼50% due to the consumption of CO 2(aq) at the cathode by reaction R6, forming carbonates (HCO 3 − /CO 3 2− ).It is also suggested that CO 2 is formed at the anode due to crossover of these carbonate anions through the membrane to the anode side and their subsequent conversion back to CO 2(aq) at the anode. 14,24inetics Governing Chemical Reactions.Reactions R1−R3b are governed by the Tafel equation: where i m is the current density of EC reaction m, i 0,m is the exchange current density.The Tafel slope A m is given by where R is the universal gas constant, T is the temperature, and a m is the anodic or cathodic kinetic fitting parameter for reaction m.The overpotential η m is given by where φ s and φ l are the potentials in the solid and liquid phases, respectively, E 0,m is the standard reduction potential, v i , m is the stoichiometric coefficient of species i in reaction m, c i is the concentration of species i, n m is the number of electrons transferred, and F is Faraday's constant.For eqs 1 and 2, the signs are positive for an anodic reaction and negative for a cathodic reaction.All EC reactions are summed over to give the total rate of production/consumption for each species i according to where a (c/a) is the specific surface area of the cathode/anode.For solution-phase reactions (reactions R4−R8), forward rate constants and equilibrium constants k j , K eq , j of each  See the Results section for a discussion on the effect of each parameter on the output of the simulation based on the results of the sensitivity analysis.
The Journal of Physical Chemistry C reaction j are supplied as input parameters, and reverse rate constants calculated according to The rate r j of each solution-phase reaction j is defined as where A−D represent activities of the species involved and a− d their stoichiometric coefficients (the activity of a species is its concentration divided by 1 M or 55 M in the case of water 4 ).
For each species, the rate of production/consumption is calculated by summing over all of the reactions it is involved in: CO 2(aq) is the reactant in reaction R1 rather than CO 2(g) since it is suggested that H (aq) + ions are needed for the reaction to occur. 7We follow the method of Weng et al. 9 and others and define within the cathode the rate of reaction R9 as where TF is the gas to liquid mass transfer coefficient, (δ TF is the thickness of the electrolyte thin-film layer coating the catalyst), w COd 2 (g) is the mass fraction of CO 2(g) , D COd 2 (aq) is the diffusion coefficient for CO 2(aq) in water, H COd 2 (aq) is Henrys constant for CO 2 in water at 298 K, p g is the gas pressure, a c is the specific surface area of the cathode, and C COd 2 (aq) is the concentration of CO 2(aq) .For CO 2(aq) and CO 2(g) , solvation (eq 8) also contributes to the rate of production/consumption.Solvation between phases is simulated only for CO 2 .
Equations Governing Transport of Species and Electrons�Mass Balance and Current Balance.Mass balance for species in the solution phase is governed by the Nernst−Planck equation: where u is the convective velocity (which is zero in this simulation), R i,total is the total rate of production/consumption of species i due to EC, solution-phase reactions, and solvation reactions: The molar flux J i of each species is defined according to the Nernst−Planck equation, describing transport due to diffusion and migration at infinite dilution: where ϵ l is the porosity of the given domain, D i is the diffusion coefficient for species i in solution, z i is its charge, and u i is its mobility according to the Nernst−Einstein relation.Increasing the porosity will allow for increased transport of liquid species but result in a decrease in the current/transport of electrons (eq 14 below) and of gases.The solid volume fraction ϵ s is calculated from the porosity as (ϵ s = 1 − ϵ l ).In the cathode, the liquid volume fraction is (ϵ l S), and the gas volume fraction is (ϵ l (1 − S)), where the saturation S = 0.64 as in Weng et al. 9 This results in the anode and membrane being 60% liquid and 40% solid and the cathode being 39% liquid, 21% gas, and 40% solid.
Electrons are not treated as a chemical species with mass; their transport is governed by the current balance in the solid phase.Current balance is governed by eqs 12 and 13: where i s and i l are the current densities in the solid and solution (liquid) phases, respectively, given by and where σ is the conductivity of the solid phase (i.e., the electrode), and φ s is the potential in the solid phase.Mass transport of gas-phase species n is governed by where ρ is the mixture density, u is the average velocity of the mixture, and w n is the mass fraction of species n, where we impose the condition ∑ n w n = 1.The mass flux j n is defined as where M is the mean molar mass , where the Stefan−Maxwell diffu- and the Knudsen diffusivity are occurring in parallel.Here, x n is the mole fraction of species n.Darcy's law governs the convective velocity according to In eqs 16 and 18, R n,total for gas-phase species is defined by eq 7 for CO 2 , and it is zero for other gaseous species.The convective velocity, μ, is defined as where p is the pressure, κ is the permeability of the medium (calculated from the Kozeny−Carman relation where the particle diameter is taken as twice the pore radius), and μ is the dynamic viscosity of the mixture.Electroneutrality is assumed in all domains of the simulation: Boundary Conditions.At the leftmost boundary (gas channel/cathode boundary), the electric potential of the solid The Journal of Physical Chemistry C phase is fixed (φ s = E applied ), and at the rightmost boundary (anode/anolyte channel boundary) it is set to ground (φ s = 0 V), corresponding to the locations of the current collectors.
At the leftmost boundary, the mole fraction of CO 2(g) is fixed corresponding to the inlet feed composition; 99% for CO 2(g) , <1% for CO (g) /H 2(g) .At this boundary, the fluxes of gaseous species can be measured, allowing us to estimate the composition of the gas channel.The flux of CO 2(g) is set to Input Parameters.The equations above are populated with parameters from Table 2 before being solved.
When solving the simulation in COMSOL, a user-controlled mesh was used with 1000 mesh elements in the anode and cathode domains combined, and 100 in the membrane domain.One mesh refinement was performed in the membrane domain.The solution was independent of an increasing number of mesh elements.The default PARADISO solver was used.Separate study steps were used to increase the membrane fixed charge to its final value and to increase the voltage from 1 to 3 V.For each step, the maximum number of iterations increased to 1000, and the relative tolerance was 1 × 10 −4 .
Sensitivity Analysis Methodology.First, a base case simulation was defined as outlined above and solved.The simulation was then solved again once for each parameter of interest, where the parameter was increased by 1% from its original value.The sensitivity, s, of the CO partial current density, i CO , to each parameter is defined as where i CO(+1%) is the CO partial current density of the simulation with the parameter increased by 1%, and i CO(base) is the CO partial current density of the base case simulation.Sensitivity analysis was also performed for a 0.1% perturbation

The Journal of Physical Chemistry C
of the input parameters, and the same sensitivity trends are reported (Figure S5 of the Supporting Information).

■ RESULTS AND DISCUSSION
Polarization Curve and Faradaic Efficiency. Figure 2a shows that using KOH as the electrolyte instead of KHCO 3 will result in a greater current density at the same potential as expected. 11,36This is due to a reduction in the Nernstian (concentration) overpotential at the anode, due to the constant supply of OH − via the anolyte from KOH. 11 OH − also has a greater ionic conductivity than HCO 3 − , resulting in less Ohmic resistance in the cell.
The FE is a measure of the utilization efficiency of the electrons in the system�it is defined as the fraction of the electrons, which go to the formation of the relevant product (CO vs H 2 ).The FE of CO is calculated as .
Figure 2b shows that at higher current densities, H 2(g) formation dominates over CO (g) formation at the cathode, as in the experiment of Larrazabal et al. 37 No significant difference is observed in the simulation for the FE of the KHCO 3 and the KOH cases at the same current densities; 11 a difference is only observed when plotting FE vs potential, since each case gives a different polarization curve.We note that the cation type and concentration are reported to affect the activity and selectivity of EC reactions, although the mechanism through which this occurs is not fully understood.One possible explanation is through the stabilization of reaction intermediates. 38With the current state of art in our knowledge of these processes, multiphysics simulations such as those presented here do not account for reaction intermediates, as discussed later, or a full set of solution-phase reactions including those involving K + .We further note that a range of different experimental Faradaic efficiencies have been reported in the literature for GDE electrolysis cells, and the EC parameters used here replicate the behaviors reported by Larrazabal et al. 37 whose setup this simulation resembles.Concentrations.Figure 3a,b shows the percentages of CO 2(g) , CO (g) , and H 2(g) in the gas channel for the 0.1 M KHCO 3 and 0.1 M KOH cases.The composition is estimated by measuring the molar flux of each species at the cathode/gas channel boundary and dividing by the total flux of all species (CO 2(g) , CO (g) , and H 2(g) ) at that boundary.In the simulation, the mole fractions are fixed at the left-hand boundary at the inlet feed composition, which is 99% CO 2(g) and <1% CO (g) / H 2(g) (i.e., the simulation assumes a sufficiently high flow rate of CO 2(g) to maintain this composition).With applied current density, the concentrations of CO (g) and H 2(g) in the gas channel increase while the concentration of CO 2(g) decreases.
For the KOH case, the CO 2(aq) being consumed at the anode in the solution phase results in diffusion out of the cathode into the membrane, and thus a lower CO 2 utilization at very low current densities is observed compared to the KHCO 3 case (see Table 3).
Figure 4a−d shows the concentration profiles of all species in all domains of the cell, for both the 0.1 M KHCO 3 and 0.1 M KOH cases, at 10 and 1000 mA/cm 2 showing the effects of applied current density and electrolyte type.CO 2(aq) is formed from CO 2(g) within the cathode, where its concentration is maintained at ∼30 mM.For the KHCO 3 case, Figure 4a,b, there is always a relatively high amount of CO 2(aq) at the anode, as HCO 3(aq) − is supplied via the anolyte, which is converted back to CO 2(aq) via reaction R4.At the anode, CO 2(aq) exits via dissolution to CO 2(g) .At higher current densities, CO 2(aq) in the membrane is also consumed by OH (aq) − , which is produced during CO 2 reduction at the cathode.For the KOH case, Figure 4c,d, we are supplying OH (aq) − via the anolyte, which consumes CO 2(aq) in the anode region via reaction R6.Thus, there is less CO 2(aq) at the anode, compared to the KHCO 3 case.This net consumption of KOH and conversion to (bi)carbonates is described by Rabinowitz and Kanan 24 and identified as a problem in electrolyzers of this type.For the KOH case, the anode contains more CO 2(aq) at high current densities as HCO at the cathode at higher potentials indicates salt (K 2 CO 3 /KHCO 3 ) formation.For the KOH case, the concentration of HCO 3(aq) − is higher at the anode at higher current densities due to the conversion of H (aq) + to HCO 3(aq) − .The concentrations of all ions are shifted in the AEM depending on their charge, as anions are attracted to and cations repelled from the membrane.

OH (aq)
− is generated at the cathode electrochemically via CO 2 reduction and is consumed at the anode electrochemically by reaction R3b.For the KOH case, the concentration of OH (aq) − is fixed at the anode boundary at 0.1 M, corresponding to a fresh supply of KOH.OH (aq) − is consumed by reaction R6, which is responsible for producing HCO 3(aq) − . 24K (aq) + is fixed at 0.1 M at the anode for both cases.K (aq) + migrates toward the cathode due to its negative charge and due to the high presence of anions there (HCO 3(aq) − , and CO 3(aq) 2− ).It is important to note that K (aq) + is not involved in any chemical reactions in this simulation.This is a simplification that could be amended by including a full set of reactions, as discussed in the Supporting Information.High concentrations of K (aq) + at the cathode at higher potentials indicate salt formation (the solubility limit is ∼8 M for K 2 CO 3 ).This occurs at approximately 750 mA/cm 2 , according to Weng et al. 11 CO 3(aq) 2− is fixed at the anode boundary at its equilibrium value for each case.CO 3(aq) 2− is produced at the cathode region primarily via reaction R7 due to the presence of OH − ions.As a For both cases, utilization (ε CO ) converges to ∼50% at higher current densities. 11The magnitudes of each process in mol/m 2 /s are included in parentheses for reference.

The Journal of Physical Chemistry C
is the case for K (aq) + , a high concentration of CO 3(aq) 2− at the cathode at higher current densities indicates salt formation.

H (aq)
+ is fixed at the anode boundary at its equilibrium value for each case.The EC oxidation of water (reaction R3a) results in a higher H (aq) + concentration at the anode.The concentration of H (aq) + is low at the cathode where OH (aq) − is produced electrochemically.Note that the H (aq) + concentration is multiplied by 10 3 so it appears higher in the figure .The concentrations of gaseous species CO (g) and H 2(g) are also shown.Although H 2(g) is produced at a higher rate at 1000 mA/cm 2 than CO (g) , the CO (g) concentration is higher since H 2(g) diffuses out of the cathode more easily, hence why we estimate the gas channel composition as in Figure 3a,b rather than paying particular attention to the gaseous concentrations within the cathode.We also note that the volume fractions in the solution phase and gas phase are different because of the

The Journal of Physical Chemistry C
porosity values used here, as discussed in the Methodology section.CO 2 Utilization.As well as current density and FE, CO 2 utilization is an important metric in determining the performance of CO 2 electrolysis cells.CO 2 utilization is a key problem for zero-gap MEA cells, as there is competition between the EC reduction of CO 2(aq) and solution-phase reactions converting CO 2(aq) to HCO 3 − /CO 3 2− . 4,13,24 Figure 5 shows a mass flux diagram for CO 2 in a typical MEA-type cell, illustrating this competition.Table 3 indicates the rates of each process ε in Figure 5; CO 2(aq) to CO (g) , CO 2(aq) consumed into solution phase (converted to HCO 3 ), and CO 2(aq) diffusing unreacted out of the cathode to the anode side.Utilization here is calculated as where r 1 refers to the reaction rate of R1 (i.e., CO 2(aq) reduced to CO), R COd 2 ,total is the total rate of CO 2 consumed, and J exit is the flux of CO 2 exiting the cathode.To find the rate of CO 2(aq) consumed, R COd 2 ,total , we sum the rates of reactions R1, R4, and R6.The conversions ε solution and ε diffusion are calculated similarly.The integrals are computed over the length of the cathode.At low current densities for the KOH case, the concentration gradient of CO 2(aq) between the cathode and anode means a high rate of diffusion of CO 2(aq) from the cathode, thus lowering both the fraction of CO 2(aq) converted to CO (ε CO ) as well as the fraction of CO 2(aq) consumed in solution in the cathode (ε solution ).We can see that for the KHCO 3 and KOH cases, the CO 2(aq) utilization ε CO converges toward ∼50% at higher current densities, while diffusion out of the cathode becomes negligible.As discussed in Weng et al. 11 and Liu et al., 39 studying the reaction stoichiometry suggests a utilization of ∼50% for these reactors (depending on the relative rates of each reaction).This upper limit is imposed by the zero-gap nature of the cells, resulting in competition between the solution-phase and EC reactions.As in the simulation of Weng et al., 11 high CO 2 utilization can be achieved at relatively low current densities of approximately 100 mA/cm 2 .This is advisable to prevent salting, flooding, and wear of materials, which all occur at higher current densities.Parameters of the multiphysics simulation can be tuned to find the optimal conditions under which to run the cell, but that is outside the scope of this paper.
Sensitivity Analysis.The results of sensitivity analysis are shown in Figure 6a,b for the 0.1 M KOH and 0.1 M KHCO 3 cases, at 10 and 1000 mA/cm 2 .We can see that, at any condition, i CO is highly sensitive to the kinetic parameters (i 0 , α) of all EC reactions.Even at high currents, i.e., in the "mass transport-controlled regime", these parameters give the biggest change in i CO (although diffusion coefficients show a higher response at higher current densities as expected).We can thus conclude that to make the greatest gains in i CO , effort should be focused toward finding better-performing catalysts with high activity toward the desired product.Indeed, this is where most work is done today in CO 2 reduction. 40The higher sensitivity to kinetic parameters of reaction R2, the hydrogen evolution reaction, is due to this reaction dominating over CO 2 reduction at high current densities.
It is important to note that the simulation is also highly sensitive to the parameters of reactions R3a and R3b, which are conventionally called the OER parameters.Whether the simulation models the OER or not and the choice of OER parameters will affect the extracted i 0 /α c for reaction R1 since many simulations exclude the anode.We also note that it is common to round these parameters, as well as kinetic parameters of solution-phase reactions, to the nearest order of magnitude from their reference source�doing so will impact the simulation results significantly as well as the extracted i 0 /α c values.
i CO also shows high sensitivity to physical parameters, e.g., the lengths and specific surface areas of the electrodes.High sensitivity of the simulation to these physical parameters indicates the need to ensure that the simulated conditions precisely match those of the experiment if one is using a multiphysics simulation to make accurate predictions about one's experiment.Further, if using a multiphysics simulation together with an experimental setup to extract EC kinetic parameters, then incorrect values will be extracted if the parameters such as lengths of simulation and experiment are not precisely matching.
Lower sensitivity is shown to rate constants for solutionphase reactions and for diffusion coefficients.However, this does not imply that these parameters are not important for the simulation, only that they do not impact i CO in particular.Other metrics of interest may be highly sensitive to these parameters.Given the scale of the uncertainty in solutionphase kinetic parameters, in particular, we recommend that further efforts are made to determine their values more precisely.
Since many parameters are difficult to determine, and simulation output is highly sensitive to many parameters, it is likely that such multiphysics simulations will not be able to make precise quantitative predictions about cell performance, e.g., predicting exact species concentrations.We recommend therefore that simulations of this kind, working with high levels of uncertainty in their input parameters, should be used to investigate trends rather than to make precise quantitative predictions, unless parameters are determined more precisely.We further note that infidelity in the simulation due to inaccurate input parameters will propagate the nature of the highly coupled reaction network in this system (Figure S6).
A sensitivity analysis of this nature can also help in determining which parameters will affect the convergence of these simulations the most if convergence issues arise due to unphysical quantities being computed.Convergence is an issue limiting the development of these models. 4,15,17,23,27By choosing a different sensitivity metric other than i CO , one could optimize cell performance by identifying which parameter most affects the metric of interest, e.g., FE or CO 2 utilization.
The Tafel equation is plotted in Figure 7 showing the current density of reaction R1 (CO 2 reduction to CO) vs overpotential for each of the i 0,1 /α c,1 pairs reported in Table 1.In the Tafel equation, i 0 (exchange current density) represents the current density at a 0 V overpotential, and α (kinetic transfer coefficient) represents the sensitivity of the reaction rate to the overpotential.We can see that each i 0,1 /α c,1 pair results in significantly different behaviors.Some disagreement in the reported values of i 0,1 /α c,1 should be expected, however, since kinetics will be affected not only by the catalyst material but also by electrolyte The Journal of Physical Chemistry C concentration, cation type, surface coverage, partial pressure of CO 2 , and morphology of the catalyst. 8,41,42The Tafel slope (or equivalently α) is also reported to change vs overpotential. 41ultiphysics simulations generally do not account for the fact that i 0 /α changes with conditions but use one number for i 0 /α in all conditions.This is another reason to expect infidelity in such simulations vs experiment.Table 1 shows that significantly different conditions are simulated in each instance.Thus, we should not expect just one i 0 /α pair for these simulations for any reaction.
Another important consideration, overlooked in multiphysics simulations, is that the CO 2 reduction reaction is modeled as a single step, when in reality it is multiple steps, each step having its own rate, and is condition dependent. 41ultiphysics simulations model this reaction as a single step using the i 0 /α of the rate-determining step.Which step is the rate-determining step, however, will change with the condition, hence using a single i 0,1 /α c,1 pair for the CO 2 reduction reaction under all conditions will result in infidelity.To accurately extract kinetic parameters i 0 /α using multiphysics simulations, one should account for the multistep nature of EC reactions, as well as condition-dependent kinetics, which is often overlooked.Attempts have been made to determine i 0,1 / α c,1 more rigorously 8,41,42 by combining multiphysics models for kinetics and transport with quantum mechanical models using density functional theory (DFT).These works account for the multiple steps in each EC reaction mechanism and how their rate constants change in different conditions.Combining a DFT model with a microkinetic and transport model will more accurately reproduce experimental results from CO 2 electrolysis cells.Singh 8 et al. use a combined multiphysics and DFT model to determine i 0,1 /α c,1 for their specific Ag catalyst in question.A multiphysics simulation alone can still be used to gain insights into cell behaviors, which are due to reaction/transport phenomena, but to build the simulation such i 0 /α needs to be chosen which accurately replicates cell behavior.

■ CONCLUSIONS
A zero-gap membrane electrode assembly cell for CO 2 reduction is simulated in 1D, considering EC and solutionphase reactions and species transport due to diffusion and migration under a zero-flow condition.The cathode is a porous Ag GDE, and the anode is a porous iridium oxide electrode supplied with aqueous 0.1 M KOH/KHCO 3 .It is shown that a significant difference in the polarization curve is observed between the 0.1 M KHCO 3 and 0.1 M KOH cases, mostly due to the constant supply of OH (aq) − at the anode for the KOH case, resulting in a more favorable anodic overpotential.At the cathode, CO formation dominates at lower current densities and is surpassed by H 2 production at current densities higher than ∼700 mA/cm 2 , which is in agreement with observed experimental trends for electrolyzers of this type.Flux analysis shows that CO 2(aq) is consumed by three competitive processes; conversion to CO, conversion to carbonates (HCO 3 − /CO 3 2− ), and diffusion out of the cathode.In agreement with previous experimental work, the CO 2 utilization metric shows that CO 2 reduction to CO is limited to ∼50% due to competition between EC reduction of CO 2(aq) and its solution-phase reaction with OH (aq) − to form the carbonates.Furthermore, if KOH is used as the electrolyte, KOH will be consumed by solution-phase reactions at the anode, forming carbonates.
Sensitivity analysis of the multiphysics simulation shows that the current density of the CO 2 reduction reaction (i CO ) is highly sensitive to the EC kinetic parameters of all EC reactions occurring: (i) CO 2 R, (ii) HER, and (iii) OER.For the CO 2 R reaction, it is further shown that the range of EC kinetic parameters i 0 and α suggested in the literature for Ag electrodes results in variance in partial current densities spanning 12 orders of magnitude.The sensitivity of i CO to many of the simulation's input parameters helps to explain the large range in the reported values of i 0 and α, which have been extracted by multiphysics simulations fitting i CO to experimental data.It is thus concluded that the veracity of multiphysics simulations is dependent on improved knowledge of the EC reactions taking place and of their chemical kinetics.However, it is noted that multiphysics simulations still offer utility if appropriately informed by experiment.Finally, the sensitivity analysis methodology presented here can also be employed to optimize cell performance, by giving the most important parameter with respect to any performance metric such as i CO , FE, energy efficiency, CO 2 utilization, etc.

Figure 1 .
Figure 1.Three domains (cathode, membrane, and anode) in the 1D simulation of the CO 2 electrolysis cell presented in this work.Below the 1D domains simulated, the MEA cell setup is illustrated.The cell consists of a cathode (a GDE of pure porous silver, i.e., pure catalyst), an anionexchange membrane, and an anode (carbon paper coated with an IrO 2 catalyst layer).At the leftmost boundary (gas channel/cathode boundary), we fix the mole fraction of CO 2(g) corresponding to inlet CO 2(g) feed.At the rightmost boundary (anode/electrolyte channel boundary), we fix the concentrations of solution-phase species H 2 O (aq) , H (aq) + , OH (aq) − , HCO 3(aq) − , CO 3(aq) 2− , and K (aq) + , corresponding to delivery of anolyte solution (either 0.1 M KHCO 3 or 0.1 M KOH).The table below the illustration shows the species and reactions in each domain, the main equations used to simulate reaction kinetics and transport, and the boundary conditions used.All equations and terms are described in the Methodology section.The equations are populated with parameters from Table2, and they are solved using COMSOL Multiphysics 6.0.Equations are described in Newman.33 u s i o n c o e ffi c i e n t D n m , k i s c a l c u l a t e d a s

r a 9 cathode
at the cathode/membrane boundary, where r 9 is the rate of reaction R9, CO 2(g) ⇄ CO 2(aq) .The fluxes of other gaseous species are set to zero at this boundary.At the anode/anolyte channel boundary, the concentrations of solution-phase species H 2 O (aq) , H (aq) + , OH (aq) − , HCO 3(aq) − , CO 3(aq) 2− , and K (aq)+ are fixed at their equilibrium concentrations corresponding to the delivery of anolyte, either 0.1 M KHCO 3 or 0.1 M KOH.These values are given in the Supporting Information.

Figure 2 .
Figure 2. (a) Total current density vs the applied cell potential and (b) Faradaic efficiencies of each product vs current density, for the 0.1 M KHCO 3 and 0.1 M KOH cases.The solid line represents the KHCO 3 case and the dashed line the KOH case.

Figure 3 .
Figure 3.Estimated percentages of species CO 2(g) , CO (g) , and H 2(g) in the gas channel vs current density, for (a) 0.1 M KHCO 3 and (b) 0.1 M KOH.The gas channel composition is estimated by measuring the fluxes of each species at the gas channel/cathode boundary.

−
profile shown in Figure 4a,b for the KHCO 3 case shows the HCO 3(aq) − concentration fixed at 0.1 M at the anode corresponding to 0.1 M KHCO 3 electrolyte used.HCO 3(aq) − is produced at the cathode via reaction R6 and subsequently is converted to CO 3(aq) 2− via reaction R7.High concentrations of CO 3(aq) 2− and of K (aq) +

Figure 4 .
Figure 4. Concentrations of all species in all domains of the cell.The solid lines represent the 0.1 M KHCO 3 case and the dashed lines the 0.1 M KOH case: (a) KHCO 3 10 mA/cm 2 , (b) KHCO 3 1000 mA/cm 2 , (c) KOH 10 mA/cm 2 , and (d) KOH 1000 mA/cm 2 .CO 2(aq) forms in the solution phase from CO 2(g) within the cathode.Concentrations of H 2 O (aq) , H (aq) + , OH (aq) − , HCO 3(aq) − , CO 3(aq) 2− , and K (aq) + are fixed at the rightmost boundary corresponding to the feed of the anolyte solution in each case.Note that the H (aq) + concentration is multiplied by 10 3 , so it appears higher in the figure.

Figure 5 .Figure 6 .
Figure 5. Flux diagram for CO 2 in the zero-gap MEA cell.The region inside the dashed border represents the processes analyzed using the 1D simulation.A key problem for this setup is competition between the EC reduction of CO 2(aq) to CO (g) and solution-phase reactions converting CO 2(aq) to HCO 3(aq) − /CO 3(aq) 2− .

Table 1 .
Range in Reported EC Kinetic Parameters i

Table 2 .
Input Parameters Provided to the 1D Simulation a

Table 3 .
Utilization of CO 2(aq) for the KHCO 3 and KOH Cases, Corresponding to the Dashed Region in Figure5 a

The Journal of Physical Chemistry C ■ AUTHOR INFORMATION Corresponding Authors Harry
Dunne − School of Physics, Trinity College Dublin, Dublin D02 PN40, Ireland; orcid.org/0000-0001-5033-631X;Email: dunneh1@tcd.ieStephen Dooley − School of Physics, Trinity College Dublin, Dublin D02 PN40, Ireland; Email: stephen.dooley@tcd.ieSchool of Physics, Trinity College Dublin, Dublin D02 PN40, Ireland Mohammad Reza Ghaani − School of Physics, Trinity College Dublin, Dublin D02 PN40, Ireland Kim McKelvey − MacDiarmid Institute for Advanced Materials and Nanotechnology and School of Chemical and Physical Sciences, Victoria University of Wellington, Wellington 6012, New Zealand; orcid.org/0000-0002-0686-6898Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jpcc.4c00690Funding This work is funded by the Sustainable Energy Authority of Ireland (Grant RDD517), by Science Foundation Ireland (Grants 12/RC/2302_P2, 12/RC/2278_2, 2/RI/9673) receiving co-funding from the European Regional Development Fund via the AMBER award and is supported by the European Union through the European Research Council's Mod-L-T project (action number 101002649).