Transient Diffusion-Reaction Model for CO2 Electrolysis

Solves the set of coupled PDEs describing mass conservation across a diffusion layer with an electrochemical reaction occurring at one side.

This model has two components explained in this order below:

• MODEL I: Transient-Diffusion Reaction (TDR) Model for a Planar Electrode
• MODEL II: The Transient-Diffusion Reaction-Kinetic (TDRK) Model for a Planar Electrode
• MODEL III: Transient-Diffusion Reaction Model for a Porous Electrode (pTDR)

The first model is based on efforts presented in the following works and adapted for the current format of interactivity:

• Gupta et al. J. Appl. Electrochem. 2006, 36, 161–172
• Resasco et al. ChemElectroChem 2018, 5, 1064–1072
• Kas et al. ChemElectroChem 2015, 2, 354–358

This has as of update 2.0 been extended to also include electrode kinetics (as MODEL II), based on work in:

• Blake et al. Electrochim. Acta 2021, 393, 138987
• Burdyny et al. ACS Sustain. Chem. Eng. 2017, 5 (5), 4031–4040

The third model is an extension of MODEL I and was further inspired by:

• Raciti et al. Nanotechnology 2018, 19, 044001
• Burdyny & Smith. Energy Environ. Sci. 2019, 12, 1442–1453
• Blake et al. Electrochim. Acta 2021, 393, 138987

It is worked on coupling MODEL III with electrode kinetics as well...

A full documentation with all sources and background of relevance and more proper descriptions is currently being written.

***MODEL I: Transient-Diffusion Reaction (TDR) Model for a Planar Electrode***

Fig. 1. Model domain with an electrochemical reaction occurring at a planar electrode.

Solving the coupled species' conservation equations across a diffusion boundary layer of length δ_DBL in contact with a planar cathode at which an electrochemical surface reaction takes place (Fig. 1). The pdepe solver is used to this end to model the behaviour of dissolved CO2, HCO3(-), CO3(2-), OH(-), and H(+). The equations that are solved for the five species' concentration profiles (both in time and space) are of the form

with c as concentration, D as diffusion coefficient, and R as source/sink term for the (bi)carbonate equilibrium reactions. The reduction reactions taking place at the cathode surface enter as Neumann boundary conditions at x = 0 and two such reactions are considered. For neutral to alkaline electrolyte the acid components of these chemical reactions may be largely ignored as can be derived from a scaling analysis of the reaction terms. Furthermore, the electrochemical reactions generate OH(-) and therefore serves to raise the pH, further motivating this simplification. This is henceforth referred to as the simplified model, though the option exists to retain the full model equations.

The model file tdr_model.m is supplied with an editor file tdr_model_editor.mlx that allows modification of pertinent parameters to study its influence on for example the concentration profiles, pH, and mass transfer limiting current density. The following parameters may be altered:

• temperature T
• pressure P
• molarity c (of CO2-saturated electrolyte)
• current density j
• Faradaic efficiency η for CO
• diffusion boundary layer thickness δ_DBL

The model file runs several subroutines (attached as separate MATLAB files) in order to determine the necessary input data for the model derived from the above input parameters. These files should be located in the same directory as the model and editor files:

• henry.m CO2 solubiltity as function of temperature and pressure
• sechenov.m CO2 solubiltity as function of ionic strength
• viscosity.m viscosity of electrolyte as function of molarity
• diffusivity.m diffusion coefficients as function of temperature
• carbonateeq.m equilibrium coefficients as function of temperature and ionic strength for (bi)carbonate equilibria
• carbonatekin.m rate coefficients as function of temperature and ionic strength for (bi)carbonate equilibria
• selfionization.m pKw value of water as function of temperature
• bjerrum.m shows the distribution of CO2, HCO3(-), CO3(2-) as function of pH given equilibrium coefficients

Examples are given of the plots returned by running the model through the editor file or data that may be derived from the results.

Fig. 2. Transient concentration and pH profiles upon imposing -5 mA/cm2 at the cathode (x = 0) for 0.1 M KHCO3.

Fig. 3. Mass transfer limiting current density in 0.1 M KHCO3 as function of diffusion boundary layer (DBL) thickness. The model result is compared with the estimate based on the formula on the right (assuming a linear concentration profile).

MODEL II: Transient-Diffusion Reaction-Kinetic Model for a Planar Electrode (TDRK)

The same approach as in MODEL I is made and calculating concentrations across the DBL still relies on the pdepe.m solver. The surface concentrations at the electrode are then used to calculate the concentration polarization alongside kinetic and Ohmic polarization. One of the key results from this model is that it can predict a polarization curve. For example, using kinetic data of a gold electrocatalyst from the work of Chen et al. (J. Am. Chem. Soc. 2012, 134(49), 19969–19972) the polarization curve could be drawn. It indicates that above -0.8 V_RHE it suffers from mass transport limitation of CO2.

Fig. 4. Polarization curves in absence (dotted) and presence (solid) of transport polarization for a gold electrocatalyst as described in the text. Note the potential scale is IR-corrected.

In order to model the electrode kinetics alongside mass transport phenomena, the kinetic parameters ought to be known for every electrochemical reaction that is desired to be included. These are the exchange current density and charge transfer coefficient, which can be found from a Tafel plot. Moreover, an appropriate Butler-Volmer expression needs to be chosen to model the kinetics in a mixed kinetic-mass transfer regime. For the time being, only H2, CO, and C2H4 are included with assumed rate expressions:

Concentrations with subscript 0 denote the bulk concentrations. One may consult the book by Oldham et al. (Electrochemical Science and Technology: Fundamentals and Applications, 1st ed.; 2012) to see how such expressions may be derived on the basis of a reaction mechanism. In the case of CO for example, its expression is based on CO2 + e(-) ---> CO2(-) being the rate-determining step, making its partial current density first-order dependent on CO2 concentration.

The general procedure of the model is shown below. Just as before, the model input requires the user to define a current density and Faradaic efficiency (for each of the products involved). In addition, the kinetic parameters should be specified for each (or set their exchange current density or Faradaic efficiency to zero to ignore). Then it proceeds through an iteration loop for the following reason:

1. current density determines reaction rate
2. rate determines electrode surface concentrations
3. concentrations determine electrode polarization (see Butler-Volmer expression)
4. polarization determines overpotential
5. overpotential determines current density

Fig. 5. General procedure in achieving self-consistency among concentration, current density, and potential.

This loop is repeated until self-consistency is achieved given by a convergence criterion C, or until the maximum number of iterations is reached. The constraint for the electrode potential E is that the sum of all partial current densities must fit the input current density j, so that

The convergence criterion is set such that the deviation in Faradaic efficiency for all products i in successive iterations (between iteration n and n - 1) is less than C, thus

To facilitate conversion at current densities close to the mass transfer limiting current density, a relaxation parameter lambda may have to be specified (between 0 and 1), defined as

The model file tdrk_model.m is supplied with an editor file tdrk_model_editor.mlx that allows modification of pertinent parameters to study its influence on for example the concentration profiles, pH, and overpotentials. The following parameters may be altered:

• temperature T
• pressure P
• molarity c (of CO2-saturated electrolyte)
• current density j
• Faradaic efficiency η for CO (initial guess)
• exchange current density j_0 and charge transfer coefficient alpha for each electrochemical reaction
• diffusion boundary layer thickness δ_DBL
• maximum number of iterations to allow
• convergence criterion C
• relaxation parameter lambda

The model file runs several subroutines (attached as separate MATLAB files) as before and a few other subroutines are introduced for incorporating electrokinetics in the model. These files should be located in the same directory as the model and editor files:

• nernst.m structure element with thermodynamics of electrochemical reactions
• butlerVolmer.m structure element with kinetic, transport, and Ohmic polarization

Both subroutines are updated during the simulation and contain the necessary kinetic data. The Butler-Volmer expressions listed above can be found in the latter and may be modified if different expressions need to be modelled. The Ohmic polarization can simply be added to a polarization curve as the product of current and high frequency resistance measured in the experiment (e.g. obtained by current interrupt method or impedance spectroscopy).

Figure 6. Kinetic, transport, and Ohmic polarization (from ionic resistance in the electrolyte) superimposed.

MODEL III: Transient-Diffusion Reaction Model for a Porous Electrode (pTDR)

Fig. 7. Model domain with an electrochemical reaction occurring in a porous electrode.

Solving the coupled species' conservation equations across a porous catalyst layer of length δ_CL and diffusion boundary layer of length δ_DBL with an electrochemical reaction taking place in the former (Fig. 4). The pdepe solver is used to this end to model the behaviour of dissolved CO2, HCO3(-), CO3(2-), OH(-), and H(+). The equations that are solved for the five species' concentration profiles (both in time and space) are of similar form as in MODEL I.

The reduction reactions take place homogeneously in the catalyst layer (0th order rate model) and enter as volumetric reaction term in the governing equations. This assumption may break down when significant gradients in CO2 concentration are found in the catalyst layer, and would necessitate a rate order showing dependence on CO2 concentration. If the associated reaction rate constant is known or specified, this may be implemented in the current model (see also the sourced work of Raciti et al.). Support is provided for the following reactions, which can be included or excluded by setting their appropriate Faradaic efficiencies.

A no-flux boundary condition is in place for all species at x = 0 with the exception of CO2 concentration (as it is supplied from the gas-diffusion layer side). This is specified as flux boundary condition with k_m as overall mass transfer coefficient. The exact value may be difficult to specify but it is likely <0.1 m/s (see also the sourced material above). At x = 0 all species are at their bulk values, dictated by the chosen electrolyte. In KOH this means that c_CO2 = 0, but in KHCO3 the bulk concentration of CO2 > 0 since HCO3(-) naturally exchanges to some degree (for as prepared 1 M KHCO3 electrolyte, about 3 mM of CO2 is naturally present according to thermodynamic data on its equilibrium coefficients).

Support is provided for neutral and alkaline electrolyte (KHCO3 or KOH) and the acid components of the chemical reactions are ignored as reasoned in MODEL I.

The model file ptdr_model.m is supplied with an editor file ptdr_model_editor.mlx that allows modification of pertinent parameters to study its influence on for example the concentration profiles, pH, and mass transfer limiting current density, or carbon efficiency. The following parameters may be altered:

• temperature T
• pressure P
• molarity c (of CO2-saturated electrolyte)
• current density j
• Faradaic efficiency η for product i
• catalyst layer porosity ε
• catalyst layer thickness δ_CL
• mass transfer coefficient k_m
• diffusion boundary layer thickness δ_DBL

These model and editor files rely on the same subroutines as MODEL I does, thus they are to be placed in the same directory as these scripts. Example below of the plots returned by running the model through the editor file.

Fig. 8. Transient concentration and pH profiles upon imposing -200 mA/cm2 to the cathode (δ_CL = 1 μm, δ_DBL = 1 μm) for 0.5 M KOH.

In addition, carbon efficiency is also calculated as either the fraction of CO2 consumed in the CL versus that entering the system at x = 0 (η_DBL) or the fraction of CO2 utilized in CO2RR compared to total consumption in the CL (η_CL). Their product is reported as total carbon efficiency. Note that this is distinct from single-pass conversion efficiency of an electrolyzer, of which this model has nothing to say about.

Fig. 9. Visual and formula representation of the two carbon efficiencies and overall carbon efficiency.

Fig. 10. Plot of carbon efficiency as function of DBL thickness for a GDE exposed to 1.0 M KHCO3 with η_CO = 0.9 showing how buffering action by the electrolyte decreases with diffusion length and lowers overall carbon efficiency.